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Galactic winds driven by cosmic-ray streaming 
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ABSTRACT 

Galactic winds are observed in many spiral galaxies with sizes from dwarfs up to the 
Milky Way, and they sometimes carry a mass in excess of that of newly formed stars 
by up to a factor of ten. Multiple driving processes of such winds have been proposed, 
including thermal pressure due to super nova- heating, UV radiation pressure on dust 
grains, or cosmic ray (CR) pressure. We here study wind formation due to CR physics 
using a numerical model that accounts for CR acceleration by supernovae, CR ther- 
malization by Coulomb and hadronic interactions, and advective CR transport. In 
addition, we introduce a novel implementation of CR streaming relative to the rest 
frame of the gas. Streaming CRs excite Alfven waves on which they scatter, thereby 
limiting the CR's effective bulk velocity. We find that CR streaming drives powerful 
and sustained winds in galaxies with virial masses M200 ^ 10 11 M . In dwarf galaxies 
(M200 ~ 1O 9 M ) the winds reach a mass loading factor of ~ 5, expel ~ 60% of the 
initial baryonic mass contained inside the halo's virial radius and suppress the star for- 
mation rate by a factor of ~ 5. In dwarfs, the winds are spherically symmetric while in 
larger galaxies the outflows transition to bi-conical morphologies that are aligned with 
the disc's angular momentum axis. We show that damping of Alfven waves excited by 
streaming CRs provides a means of heating the outflows to temperatures that scale 
with the square of the escape speed, kT oc v% sc . In larger haloes (M200 ^ lO n M ), 
CR streaming is able to drive fountain flows that excite turbulence, providing another 
means of heating the halo gas. For halo masses M200 ^ 1O 1O M , we predict an ob- 
servable level of Ha and X-ray emission from the heated halo gas. We conclude that 
CR-driven winds should be crucial in suppressing and regulating the first epoch of 
galaxy formation, expelling a large fraction of baryons, and - by extension - aid in 
shaping the faint end of the galaxy luminosity function. They should then also be 
responsible for much of the metal enrichment of the intergalactic medium. 

Key words: cosmic rays - galaxies: formation - galaxies: evolution - galaxies: star- 
burst - galaxies: dwarf - intergalactic medium 



1 INTRODUCTION 

Large-scale galactic outflows constitute an important pro- 
cess in the evolution of galaxies. They influence the gas 
content and therefore the star formation rate in galaxies. 
Moreover, outflows from galaxies also affect the surround- 
ing medium by de positing mass, metals and energy (see 
IVeilleux et al.ll2005l . for a recent review). 

The earliest works on galactic outflows were moti- 
vated by the observations of the starburst galaxy M82 
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(iLvnds fe Sandaedll96S iBurbidge et al.lll964T ). and the de- 
tection of broad optical emission lines in s ome elliptica l 
galaxies (|Osterbrockll 19601 ) that implied winds (jBurkdl 19681 ) . 
It was suggested that galactic winds cause ell ipticals to shed 
their interstellar mediu m (ISM; iMathews & Baker! Il97ll : 
I Johnson fc Axfordlll97lh . Galactic outflows have since been 
frequently observed in both local and high redshift galax- 
ies. Mult i- wavelength studies for local galaxies have demon- 
st ra ted the ubiqui ty of the phenomenon of outflows (e.g., 
iHeesen et al.l 1 20091 ) , while spectroscopic studies of Lyman 
Break Galaxie s at z ~ 3 have also r evealed the presenc e 
of winds fe.g.. lAdelberger et al.ll2003l : IShaplev et alJl2003h . 
Large mass-loading factors (defined as mass loss rate di- 
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vided by star formation rate, SFR) with values of a few 
are estimated for local galaxies and massive star- forming 
regions at z ~ 2 ! - 3 (lMartinlfl999l : iHeckman et alJlJOOOt 
ISato et aLllJKJI : ISteidel et al.ll20ld : ICoil et al.ll2oTlh . The 
use of classic ionization diagnostic diagrams allows to cleanly 
separate winds that are d riven by either starbursts or ac - 
tive galactic nuclei fAGN: ISharp &; Bland-Hawthornll2010l ): 
while AGN photoionization always dominates in the wind 
filaments, this is not the case in starburst galaxies where 
shock ionization dominates, suggesting that the latter is an 
impulsive event. Moreover, after the onset of star formation, 
ther e is a substantial delay before a starburst wind devel- 
ops ([Sharp fc Bland-HawthornlbOlOl ) . This may suggest that 
different physical processes, exhibiting different timescales, 
are responsible for driving winds in starbursts and AGNs, 
and reinforces the importance of studying different physical 
mechanisms that could drive winds. 

Theoretically, it has proven to be challenging to launch 
those powerful winds with mass-loading factors of a few 
from the galactic disc. The first models of outflows invoke 
the heating of the ISM by repeated supernovae ([Larson! 
Il974l : iDekel fc Silkl Il98d ). in which thermal pressure of 
the hot gas drives the o utflow. Early simulation studies 
([Mac Low & Ferrarall 19991 ) that placed supernova explosions 
in simple models for dwarf galaxies predicted comparatively 
little mass loss. Recently, simulations have begun to ap- 
proach the necessary resolution to explicitly simulate the for- 
mation of a multiphase medium with a dense, cold molecular 
phase that harbours star formation in full three-dimensional 
models. The collective action of stellar winds and supernovae 
drives hot chimne ys and superbubbles that start to launch 
galactic outflows (jCeverino fe Klvpinl | 2009[). which may b e 
important for low-redshift galaxies ((Hopkins et al.l l201lh . 
However, these thermal pressure-driven winds suffer from 
their tendency to destroy the cool gas through evaporation 
which limits t he mass- loading fac t or and calls for alter native 
mechanisms (JNath fe Silkl 120091 : [Murray et al] 1201 lh . Ob- 
served correlations between wind speed and rotation speed 
of galaxies, as well as the SFR, may be easier to explain with 
radiat ion pressure driving on dust grains embedded in out - 
flows (lMartinll2005l : lMurrav et~alll2005l : iNath fe Silkl l2009h . 
This appears to be one promising mechanis m of explaining 
outflo ws in high-redshift massive starbursts ((Hopkins et al.l 

[mil). 

CRs can also drive a large-scale outflow if the cou- 
pling between h igh ene r gy pa r ticles and thermal gas is 
strong enough (llpavichl Il975l; iBreitschwerdt et al.l Il99ll ; 



IZirakashvili et al] ll996J : IPtuskin et al.l Il997h . This is ex- 
pected given the near equipartition of energy in the Milky 
Way disc between magnetic field s , CRs. and the therma l 
pressure (jZweibel fe Heilesl Il997l : iBeckl l200ll : ICoxl l2005h . 
Fast streaming CRs along the magnetic field excite Alfven 
wave s through the 'streaming instability' ([Kulsrud & Pearcd 
|l969|). Scattering off this wave field limits the CRs' bulk 
speed. These waves are then damped, effectively transfer- 
ring CR energy and momentum to the thermal gas; hence 
CRs exert a pressure on the thermal gas by means of scat- 
teri ng off Al f ven w aves . 

llpavichl (|1975r ) studied this process assuming a spheri- 
cal geometry and that the waves were c ompletely damped 
away. Later, Bre itschwerdt et al.l (|l99lh considered a disc 
geometry, and calculated the effect of both small and large 



damping. They found solutions of the outflow equations 
with realistic (~ 1 Mq yr _1 ) mass loss rates from Milky 
Way-type galaxies. Such winds can explain the small gra- 
dient in 7-ray emission as a fu nction of galacto- centric ra- 
dius ([Breitschwerdt et alj|2002l ). Similarly, an outflow from 
the Milky Way driven by CRs explains the observed syn- 
chrotron emission as well as the diffuse soft X-ray emission 
towards the Galact ic bulge region much better than a stati c 
atmosphere model ([Everett et al.ll2008l : [Everett et al.ll2010h . 
This match is particularly striking in the hard X-ray band 
because of the CR Alfven- wave heating that partially coun- 
teracts the adiabatic cooling of the thermal gas and makes 
a strong case in favour of CR physics being a necessary in- 
gredient in understanding galactic winds. CR-driven out- 
flows, in which the CR fluid provides an additional source 
of pressure on thermal gas, may eject substantial amounts 
of gas from spherically symmetric galaxies with a mass out- 
flow rate per unit S FR of order 0.2 — 0.5 for massive galaxies 
([Samui et al.ll2010l ). Hence, CR pressure in starburst galax- 
ies would provide a negative fee dback to star formati on and 
eventually limit the luminosity ([Socrates et alj|2008h . 

These general ideas and analytical arguments have 
been supported by detailed numerical simulations employ- 
ing an implementation of CR physics in the Gad get code 
(ljubekas et al.ll2008l : IWadepuhl fc SpringdfeOllh . Accord- 
ing to these simulations, CRs that are accelerated by super- 
nova (SN) shocks can exert enough pressure to significantly 
influence star formation, particularly in low mass galaxies so 
that the observed faint-end of the satellite luminosity func- 
tion in the Milky Way can be successfully reproduced. These 
simulations assumed the CR fluid to be perfectly coupled to 
the thermal gas and did not allow for CR streaming in the 
rest frame of the gas. In this paper, we improve upon these 
earlier studies and investigate the effect of CR streaming 
in hydrodynamical simulations of galaxy formation. We will 
show that CR streaming is potentially the dominant energy 
transport process powering outflows from dwarf galaxies. In 
Section [2J we introduce the concept of CR streaming and 
work out a spherically symmetric wind model driven by CR 
streaming in Section [3l In Section [H we introduce our simu- 
lation setting. In Section [5] we study how CR streaming can 
drive outflows and assess their specific properties including 
mass loss rates of galaxy haloes as well as the CR streaming- 
induced heating of the halo gas. We conclude in Section [6] 
Supporting material on the numerical implementation, code 
and resolution tests are shown in the Appendices [A] and [B] 



2 COSMIC RAY STREAMING 

CR proton populations have several properties which could 
help galaxies to launch strong winds. First of all, supernovae, 
being one of the main energy sources of the ISM, are believed 
to convert a significant fraction (10% — 60%) of their released 
kinetic energy into CRs via diffus i ve shock accelerat ion in 
SN remnants dKang fc Jonedl2005l : iHelder et aljl2009h . The 
CR population forms a very light fluid in the galaxy, with 
a significant pressure and the tendency to buoyantly escape 
more rapidly than any other of the ISM components from 
the galactic disc. Since the CR fluid is coupled via magnetic 
fields to the thermal ISM components, the CR pressure gra- 
dient helps to lift thermal gas out of galactic discs. 
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Second, the energy loss time-scales of trans- and ultra- 
relativistic protons due to Coulomb and hadronic interac- 
tions with the thermal ISM are long. The CR loss times are 
longer than typical residence times of these particles in the 
denser galactic discs, where the target densities are highest 
(and therefore loss times are shortest), and they are sub- 
stantially longer than the cooling time of the thermal gas. 

Third, nearly all the energy lost by the CR population, 
be it via particle-particle interactions, via volume work, or 
via collective plasma effects (as discussed below) is deliv- 
ered to the thermal plasma. Dissipating CR energy into the 
thermal plasma increases the total pressure of the compos- 
ite fluid due to the harder equation of state of the thermal 
plasma. It can be easily verified that the increase of ther- 
mal pressure is larger than the decrease in CR pressure by 
evaluating the amount of energy transferred per volume el- 
ement, AE/V = P/(7 - 1) = P cr /(7cr - 1). The pressure 
can double in the limiting case of a fluid initially composed 
of an ultrarelativistic CR population with j C r = 4/3 which 
dissipated all its energy. More importantly, any energy in- 
jected from CRs into the thermal plasma above the galactic 
disc is protected against the severe radiative losses any ther- 
mal gas element would have suffered within the dense and 
high-pressure environment within the galactic disc. 

Thus, CRs could play an important role as a relatively 
loss-less vehicle for the transfer of SN energy into the wind 
regions. They can energize the wind with a higher efficiency 
than the thermal gas, which radiatively looses a large frac- 
tion of the energy injected into it within the galactic discs. 
The mobility of CRs to travel rapidly along already opened 
magnetic field lines and their ability t o open ISM m agnetic 
field lines via the Parker instability (| Parker! ll966T ) render 
them ideal for this function. 

The dominant CR transport process along such opened 
field lines should be streaming as we will explain now. This 
process appears naturally provided there exists a gradient 
of the CR number density along a field line. Since the indi- 
vidual CR particles travel with nearly the speed of light on 
spiral trajectories guided by the magnetic field direction, a 
gradient of CRs at a location leads immediately to a bulk 
motion of particles from the dense to the less dense region. 
This bulk motion is mainly limited by resonant scattering of 
CR particles on plasma waves. The scattering isotropizes the 
CRs' pitch angles, and thereby reduces the CR bulk speed. 

The amplitude of plasma waves which are able to res- 
onate with CRs of a certain energy therefore determines the 
rate at which the streaming CR population is re-isotropized, 
and thus the limiting bulk speed. The plasma waves are cer- 
tainly seeded by the unavoidable turbulent cascade of kinetic 
energy from large to small spatial scales in the violent ISM. 
They are further amplified by plasma instabilities, which are 
powered by the free energy of the anisotropic CR popula- 
tion. Thus the CRs, while streaming through the plasma, 
excite the waves which limit this streaming. The growth of 
the waves itself is limited by plasma physical damping pro- 
cesses. The energy of the CRs transferred to the waves is 
finally thermalized and heats the plasma. 

The delicate question is at which level does this inter- 
play of CR streaming due to a CR gradient, plasma- wave 
excitation due to the CR streaming, and wave-damping pro- 
cesses saturate? What is the resulting effective bulk velocity 
achieved by the CRs? 



A unique answer to these questions is difficult. In 
a low-beta plasma, where the Alfven speed exceeds the 
sound speed by a large factor, most works on this topic 
seem to agree that the effective streaming speed should 
be of the order of the Alfven velocity (e.g., IWentzell 1968 ; 
Kulsrud fe Pearcdll969l : lKulsrud fe Cesarskvl!T97ll ; iKulsrudl 
2005T ). In a high-beta plasma, where the Alfven speed is 
strongly subsonic, the Alfven speed can not be the limiting 
velocity for CRs. This should become clear by the following 
Gedankenexperiment. If magnetic fields become weaker and 
weaker in an otherwise unaltered plasma, the Alfven velocity 
would approach zero and any CR streaming, which would be 
limited to this speed, would vanish. However, weaker mag- 
netic fields imply weaker coupling of the CRs to the thermal 
plasma and therefore intuitively one would expect a larger, 
or at least constant transport velocity that is independent 
of the magnetic field strength in this limit. 

Consequently, a number of authors have worked 
on the problem of CR streaming in h igh-beta plasmas 
(iHolman et al.1 Il979l: lAchterberd 1 l98ll : iFelice fe Kulsrudl 
200 lj; lYan fe Lazarianll2008rh and a summary of these works 
can be found in lEnBlin et al.l (|201lh . The consensus among 
most of these authors seems to be that the limiting velocity 
is of the order of the sound speed in this situation. It can 
be several times this speed, or less, depending on details of 
the pre-existing wave level, thermal and CR energy density, 
magnetic topology and so forth. 

Here, we adopt the pragmatic approach of lEnfilin et al.l 

(|201lh to parametrize our lack of knowledge on the precise 
dependence of the magnitude of the CR streaming speed t; s t 
on the plasma parameters and assume that it is proportional 
to the local sound speed c s , 



fst 



Ac s 



(1) 



with a proportionality constant A ^ 1. 

An important aspect of this work is the energy transfer 
of the CRs to the thermal plasma. This energy transfer is 
able to repower the thermal wind continuously, in particu- 
lar at large heights above the galactic disc. The CR stream 
has a significantly higher velocity than the thermal wind, 
since the CRs' bulk speed with respect to the thermal gas 
is v B t ~ c s . Thus the dominant energy flow of a galactic 
wind region might be carried by the mobile CR fluid, and 
not by the much slower thermal gas. The transfer of en- 
ergy from CRs to the gas is mediated by plasma waves. To 
the CR population, the process of wave excitation resembles 
adiabatic work done during the expansion of the CR fluid 
with respect to the rest frame of magnetic irregularities of 
interacting plasma waves. Those are generated by stream- 
ing CRs and provide efficient scattering partners which are 
pushed away from the regions with high CR pressure. If the 
plasma wave motion would be reversed, they would com- 
press back the expanded CRs and thereby return exactly 
the transferred energy. However, since the CRs are driving 
these waves, this usually does not happen under the circum- 
stances we are interested in. Only the wave damping itself 
is a dissipative process, which thermalizes this energy and 
therefore re-powers the plasma. 

For these reasons, the CR energy loss term in the wind 
equations due to CR streaming formally resembles an adia- 
batic loss term, which would be reversible if the streaming 
velocity field converged. Just because the relevant stream- 
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ing velocity is always diverging away from the region with 
high CR pressure, the corresponding energy change of the 
CR population is always negative and the thermal plasma 
is always heated by CR streaming. 



CRs at infinity, and is negligible for an adiabatic flow. We 
are therefore left with the following expression, 



vl + 3cg )b + ■ 



'(/cr-tcr,b 
qPcr,b 



2A$. 



(7) 



3 ANALYTICAL ESTIMATES 

We begin with analytical estimates of the outflows excited by 
the heating from CR streaming. For simplicity, we assume a 
steady state situation and assume that CRs first diffuse out 
to a height above the disc, comparable to the scale height. 
Assuming spherical symmetry, we can write the constant 
gas mass flux per unit solid angle, q = pvr 2 , where p is 
the gas mass density, v is the gas velocity, and r is the 
radial distance. The mass density of CRs, p cr , is fixed by 
the constant CR flux per unit solid angle in the steady state: 
<?cr = pcrVcrr 2 , where v cv — v-\-c s is the CR speed (for A ~ 1 
in equation QJ). 

We will first calculate the terminal wind speed using 
the Bernoulli theorem, assuming that there are streamlines 
along which gas can travel from the disc to a large distance. 
In the steady state and for spherical symmetry, the energy 
equation for a compressible fluid is given by 



Or 



pvr 



v 2 P 
2 p 



F + V-Q, (2) 



where F and Q represent the external force and energy flux 
respectively, e is the specific internal energy (thermodynamic 
energy per unit mass), and P is the pressure. Here we have 
a two-component fluid composed of gas and CRs. Using an 
adiabatic index of 7 = 5/3 for the gas, and writing c 2 = 
jP/p, we have the following gas energy equation, 



r z or 



A, 



cr — heating 5 



(3) 



where 4> is the gravitational potential and A cr _ heating is the 
heating term due to damping of CR-excited Alfven waves. 
We neglect radiative cooling for the tenuous gas in the wind 
in this analytical treatment. Since p cr <C p, we can safely 
neglect the bulk kinetic energy of CRs as well as their grav- 
itational attraction. Thus, we have (using e = P cr /3 as an 
approximation for ultra- high-energy particles), 



r 2 dr \ p c 



y£L__ / a x cr 1 _ _ a 

o ^ I ^ I — -^cr — heating- 



(4) 



Here, the negative sign of A cr _heatin g indicates the loss of 
CR energy due to wave excitation. Adding these two equa- 
tions and integrating, we get the following equation for total 
energy, 



^ +3 f + *) + 4^=C, 



(5) 



where C is a constant. Equating the values at the base and 
the end of a streamline, we have 



00 1 / q s,oo 



-f^ 3 if + 



Qqcv-L cr,c 



vl c S)b 4g cr P cr , b 



+$00 = ^+3^+ 



2 2 qp C r,b 



+$b 

(6) 

The sum of the terms inside the parenthesis is 3/2 times 
the square of the effective combined sound speed of gas and 



where A$ = <£>oo - <£b- 

We will use the gravitational potential due to a baryonic 
component and the dark ma tter h alo, given by the Navarro- 
Frenk- White (NFW) profile (|l996h . Since $00 = 0, we obtain 



Vb 



GM 2i 



ln(l + rb/r s ) 



ln(l + c) - c/(l + c) 






(8) 



Here, Mbar is the total baryonic mass and M200 is the halo 
mass of a sphere enclosing a mean density that is 200 times 
the critical density of the universe, r s is the characteristic 
radius of the NFW profile, in which the density profile is 
given by p oc (r/r s ) _1 (l + r/r s )~ 2 . Assuming the initial gas 
speed ^b ~ c s ,b, the gas sound speed at the base, we have 
an expression for the wind speed at large distance which 
depends on the CR parameters and the gas sound speed at 
the base. The condition for an unbound flow or the existence 
of the wind is given by, 



4cg b + 8gcrPcr - b >2A$. 

<?Pcr,b 



(9) 



Multiplying equation (0 with the mass density, we arrive at 
the energy budget of the problem, 



£wind 



PVc 



2 4pq CY P cr , h 

: £gain - £loss = 2pC sh H p/\<P , 

pcr,bq 



(10) 

which essentially states that the kinetic energy of the wind 
is the difference between the total energy gain and the loss 
of energy against gravity. In other words, we have, 



£wind 



qA$ 



2< b + ^f 



(11) 



Now, using the value of CR mass flux at the base q c 
Pcr^b (^b + c s ,b) ~ 2p cr c s , h rl, we have 



4c 2 b + I^IV^_2A$, 



£wind 



1 



q 

qA$ 



2qc 2 h + 8c s , b rgP cr , b 



(12) 
(13) 



In the next step, we want to assess how the wind speed 
and mass- loading factor scale with the halo mass and how 
this compares to the simulations that we will introduce in 
Section U To this end, we have to specify the wind parame- 
ters which should reflect realistic ISM values and also match 
the late-time behaviour of our simulations with an approxi- 
mate steady state. We use c s ,b = 10 km s _1 , corresponding 
to an ISM temperature of 10 4 K and assume P cr ,b = 10 -12 
erg cm -3 . We adopt a mass flux (per unit solid angle) of 
q = 0.01/(47r) M0 yr _1 , which is typical in our simulations 
that show an almost spherical outflow. 

The scaling of the base height rb exhibits two 
regimes separated by a characteristic circular velocity of 



^c, crit 



120 km s" 1 . At 



^c, crit? 



the vertical scale 



height increases roughly as rb oc v c according to two- 
dimensional fits of edge-on exponential discs in the if s -band 
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Figure 1. Wind speed as a function of halo mass (equation (|12j0 
for four models where we vary the input parameters, namely the 
halo concentration c, the constant of proportionality connecting 
the base height r^ and disc scale rd, and the maximum base 
height rb,max- The blue dotted curve shows the maximum ro- 
tation speed as a function of halo mass. The triangles are the 
results of the simulation: v w i n d — (45, 135, 230) km s _1 for 
M200 = (10 9 , 10 10 , 10 11 )h~ 1 M©, while the wind in the last 
case (grey triangle) was only sustained for t < 1.5 h~ x Gyr. Filled 
circles show the circular velocity measured in the simulations at 
the virial radius. 



(|Dalcanton et al.l 12004) I 1 ! At v c > f c ,crit 3 the vertical scale 
height appears to be independent of v c . 

This behaviour can be understood theoretically by the 
following line of arguments. In shallow gravitational poten- 
tials, as found in dwarfs, the vertical scale height should scale 
with the system size. Hence, we assume the base height rb 
to scale linearly with the disc size ra (which is in turn a 
funct ion of halo mass, according to the model bv lMo et al.l 
1 19981 ) . When the gravitational potential of the halo is deep 
enough, the vertical scale height is set by the effective equa- 
tion of state of the ISM. For an effective equation of state 
of P e ff oc p 2 , the solution to the hydrostatic equation yields 
a vertical scale h eight that is i ndependent of surface mass 
density (see, e.g. , ISpringell l2000). It turns out, that the sub- 
resolution model employed by our simulations has such a 
stiff equation of state just above the star fo rmation threshold 
(see Fig. 2 of ISpringel &: Hernquistll2003h . Hence, the ver- 
tical scale height should be set by radiative cooling (which 
breaks the scale- invar iance of gravity seen in smaller sys- 
tems) in haloes heavy enough to support a thin disc, i.e., for 
M > 10 11 M© corresponding to v c > v C: C rit- This justifies a 
constant vertical scale height for these large haloes. For our 
main models, we adopt the following simple prescription for 
the base height, rb = 0.5 rd for rb < Tb, max- We note that 
these considerations are reflected by our simulations where 
we find identical trends of rb with halo mass. 

In Fig. [1] we show the wind speed at a large distance 
as a function of halo mass. This is compared with the ro- 

1 Here we assume that the base height of the wind scales with 
the vertical scale height of the disc. 
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Figure 2. Mass loading factor 77 = M/SFR (equations (|15|) and 
(1161) ) as a function of halo mass for two values of the gas fractions, 
/gas =0.1 and 0.01 (solid and dashed). The two different models 
both have a halo concentration c = 12, as in Fig. [T] The analytic 
models compare well with our simulated values for the "peak mass 
loading" (?7peak = M max /SFR max , filled circles). 



tat ion speed of th e disc imp l ied b y the halo mass, adopt- 
ing the model of iMo et all (J1998I ). We note that in this 



model, the disc mass (which we equate with Mb ar ) is a fac- 
tor ~ 0.05 of the halo mass. In the low-mass regime, the 
wind speed increases with halo mass (in accordance with 
the simulations) up to wind velocities of tVind = ^00 ~ 
(200 — 350) km s _1 , depending on the adopted parameters. 
At halo masses around M200 ~ 1O 11 M , the CR-driven 
wind is not any more powerful enough to overcome the in- 
creasing potential difference of the halo and ceases. It is 
informative to contrast this behaviour of our main models 
(solid lines) with two limiting cases where we either adopt 
a constant rb (dash-dotted) or do not impose a cutoff value 
rb,max to the rb — Td scaling (dashed). 

The behaviour of our main models is in reassuring 
agreement with our simulations where the wind speed in- 
creases with halo mass up to a critical halo mass of M200 ~ 
10 11 M©. At this mass scale, we see an onset of a wind that 
stalls after < 1.5 /i -1 Gyr, coinciding with the time of the 
formation of the thin disc. Here the increasing stellar mass 
deepens the central potential to the point where the CR- 
driving is not any more powerful enough to support a wind; 
instead we end up with a fountain flow. The normalization of 
the wind speed and the transition mass varies with changes 
in the base height, rb, and the halo concentration, i.e., the 
effective potential difference from the starting point of the 
wind to infinity. Interestingly, in the model with the high 
concentration c — 21 in Fig. [1] the wind speed starts to 
flatten towards lower halo masses. This is because for these 
small haloes the (constant) sound speed starts to dominate 
over the terms accounting for the loss of energy against grav- 
ity and the energy gain due to CR pressure in equation (p"2|) . 

To estimate the mass loading of the wind, we first es- 
timate the disc surface density Ed = Md/(7rrd), where the 
disc mass is Md = 0.05M200- The gas surface density can be 
written as E gas = / gas Sd, where / gas is the gas mass frac- 
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tion. To determine the SFR per unit area, XIs fr, we use the 
Schmidt- Kennicutt relation (|Kennicuttl ll998) , 



^SFR 



2.5 x 10" 



M© yr- 1 kpc -2 V M Q y r_1 P c ~ 2 

Hence the SFR of the disc is obtained as 



(14) 



SFR = TrraSsFR = 



2 

7rr d 



2.5 x 10" 



Jgas^d 



M, 



yr A pc" 



(15) 

We can estimate the value of the mass flux per unit solid 
angle, q, for escaping winds, by requiring that t; w ind = ^esc = 
3v c in the wind equation (equation (|12|> ) and by inverting it 
to get M = 4nq (for a spherical outflow): 

2 6 7TC s ,brbPcr,b 



M 



9^2 + 2A0 - 



4c? ' 

S,0 



(16) 



This is reasonable since the es cape speed for a NFW halo is 
< 3v c (|Sharma fc Nathll2Q12h . We can then divide this by 
the corresponding SFR of equation ([15]) to obtain the mass 
loading factor 77 = M/SFR. 

In Fig. [21 we compare the estimated mass loading fac- 
tor with the corresponding values obtained in the simulation. 
We find that our simple theoretical estimate nicely describes 
the mass loading observed in simulations, which scales as 



77 oc M : 



-2/3 



Hence we conclude that the mass loading 



in CR-driven winds scales with the halo mass in exactly the 
way that is needed to explain the low-mass end of luminos- 
ity functio n according to sem i- analytical models of galaxy 
formation (IBower et al. 20 1l|) and ph enomenological simu- 
lation models (|Puchwein &; Springelll2012r ) . CR driven winds 
are thus also very promising candidates to help explaining 
the "missing satellites problem" in the Milky Way. 



4 NUMERICAL SIMULATIONS 

4.1 Simulated physics 

We now turn to numerical simulations of galaxy formation 
that include CR physics. Our simulations were carried out 
with an updated and extended version of the distributed - 
memory parallel TreeSPH code Gadget-2 (|Springelll2005h . 
Gravitational forces were computed using a tree algorithm. 
Hydrodynamic forces are computed with a variant of the 
smoothed particle hydrodynamics (SPH) algorithm that 
conserves energy and e ntropy where appropriate, i .e., out- 
side of shocked regions (jSpringel & Hernquistll2002r ). 

Radiative cooling was computed assuming an optically 
thin gas of primordial composition (mass- fraction of X = 
0.76 for hydrogen and 1 — X = 0.24 f or helium) i n coll i- 
sional ionization equilibrium, following iKatz et al.l ([1996). 
We accou nt for star formation a c cordin g to a sub-resolution 
model by ISpringel fc Hernquistl ([20031 ) which assumes that 
star forming regions establishes a self-regulated regime, de- 
scribed by an effective equation of state. The self-regulation 
arises due to the interplay of cold, dense clouds that contin- 
uously form stars and supernova feedback that evaporates 
these clouds, creating a hot ambient medium. The rate of 
star formation is adjusted such as to reproduce the observed 
Schmidt-Kennicut law. As far as the important parameters 
of this model are concerned, we choose £q — 3.0 Gyr, thereby 
setting the characteristic time-scale of star formation, and 



a density threshold of star formation of p t h — 0.1cm -3 (in 
units of hy drogen atoms per unit volu me) . This choice is mo- 
tivated bv lDubois &; Tevssierl ([2008) who studied wind for- 
mation using similar initial conditions and a similar descrip- 
tion of star formation. Their work differs from ours in that 
they considered superbubbles as the mechanism that drives 
outflows and they observed the strongest winds for the star 
formation time-scale stated above. Since CRs as well as su- 
perbubbles are powered by supernova explosions, adopting 
their time-scale should also give the most favourable results 
in our case. 

We model CR physics with a formalism that follows 
the most important CR injection and loss processes self- 
consistently while acc ounting for the CR p r essure in the 
equat ions of motion (|Pfrommer et al.l 120061 ; lEnfilin et al.l 
2007; ljubelgas et al.l l2008). In our methodology, the non- 
thermal CR population of each gaseous fluid element is ap- 
proximated by a simple power law spectrum in particle mo- 
mentum, characterized by an amplitude, a low-momentum 
cut-off, and a fixed slope a = 2.5. Adiabatic CR transport 
processes such as compression and rarefaction, and a num- 
ber of physical source and sink terms, which modify the CR 
pressure of each particle, are included. As the most impor- 
tant source of galactic CRs, we consider CR acceleration by 
supernova remnants with a spectral index for the freshly ac- 
celerated CR population of q;sn = 2.4 and an acceleration 
efficiency, i.e., the fraction of supernova feedback energy di- 
rected into the CR population, ranging in between £sn = 0.1 
and 0.3 (where the latter value is our standard choice). As 
sinks of CRs, we consider thermalization by Coulomb inter- 
actions and catastrophic losses by hadronic interactions. 

In addition to the advective CR transport, we have im- 
plemented CR streaming relative to the rest frame of the gas. 
The details of the numerical implementation in Gadget- 2 
are described in Appendix [A] We set the streaming speed 
equal to the local sound speed (i.e., A = 1 in equation QJ). 

4.2 Initial conditions and settings 

In our simulations, we consider disc galaxy formation in iso- 
lated dark matter haloe s, according to the standard picture 
([Fall &: Efstathioulll980h . Initially, the dark matter potential 
well is filled with a hydrostatic, slowly rotating, baryonic gas 
in equilibrium with the dark matter. Right after the start of 
the simulation, radiative cooling sets in, removing thermal 
pressure support of the gas atmosphere which will then col- 
lapse, subject to gravity, to form a rotationally supported 
galactic disc in the centre of our halo. The gas in the disc 
becomes colder and denser until it starts to form stars at 
a rate given by our sub-resolution model. This picture of 
galaxy formation in isolation is clearly oversimplified, since 
galaxy formation is a highly hierarchical process, shaped 
by frequent mergers of different dark matter haloes. How- 
ever, isolated haloes enable us to study the impact of CRs 
in a comparatively clean and well- controlled environment, 
so that our simulations should be able to reveal the signif- 
icance of CR streaming for driving outflows and to make 
some qualitative predictions on the size of the effects. 

Both dark matter and baryonic gas are set up in such a 
way that they follow the same NFW density profile. Further- 
more, the two components are initially in equilibrium, i.e., 
the gas atmosphere will be stable in the absence of radia- 
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17 


(2/1) x le5 
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64 
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64 
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75 
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4.3 x 10 5 
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0.226 


0.153 


0.01 


64 



Table 1. Simulation parameters adopted in this study. Here, M200 and R200 denote the virial mass and radius, raDM and ra ga s denote 
the particle masses used to sample dark matter and baryonic gas, e^ r ^ and efr| v are the gravitational softening lengths for the two 
components, e regulates the size of the streaming time step (see equation (IA24[0 and NGB is the number of SPH smoothing neighbours 
used in this study. 



tive cooling when evolved over a Hubble time. In addition, 
the halo carries an initial amount of angular momentum de- 
scribed by the spin parameter A sp i n = 0.05, and with a radial 
distributio n in agreement with r esults from full cosmological 
simulations iBullock et al.l (|200lh . We adopt a concentration 
parameter of the NFW profile of c — 12 in all our runs. 
The initial mass fraction of baryons is kept fixed in all cases 
at its universal val ue, Ob/^m — 0.133, to ease compari- 
son to the results of ljubelgas et al.l (| 20081 ). Our simulations 
represent the baryons with 10 5 SPH-particles and the dark 
matter with twice as many particles, which produces nu- 
merically converged results (see Appendix IB2[) . Further, we 
use 64 SPH smoothing neighbours in all our simulations. 
An overview of the most important numerical parameters 
adopted in our simulations is given in Table [1] 



5 RESULTS 

5.1 Launching of galactic winds 

Since outflows cannot occur in simulations without CR feed- 
back (or radiation pressure) , at least for the specific star for- 
mation model that we use, we will in the following compare 
simulations with only advective CR transport to simulations 
with both, advective and streaming transport of CRs. Since 
the advect ion- only c ase was shown not to drive outflows 
([jubelgas et al.ll2008i ). this numerical setup enables us to 
assess the question whether CR streaming is able to launch 
galactic winds. In Figure 03 we show maps of the baryonic 
gas density in haloes of total mass 10 9 h~ x M© (top pan- 
els) and 10 10 /i -1 M© (bottom panels), seen in an edge-on 
slice through the midplane of the galaxy that forms in the 
centre. Those maps are supplemented by velocity vectors 
of the gas and are shown at t = 6/i _1 Gyr after the start 
of the simulation when the initial burst of star formation 
has ended and the galaxies have settled approximately into 
steady state. For reference, the virial radii of the two haloes 
are R200 ~ 17 /i -1 kpc and R200 ~ 35/i _1 kpc, respectively. 
When comparing the simulations with CR streaming 
and advection (left) to the advection-only case (right), it 
becomes clear that the additional mobility of the stream- 
ing CRs enables the driving of a powerful outflow. This 
manifests itself in a shallower density profile in the CR- 
streaming case, suggestive of a substantial mass transport 
beyond the virial radius. The outward pointing gas veloc- 
ities with maximal values of ^ ma x > tw demonstrate the 
presence of a wind driven by streaming CRs that extends 
well beyond 3i?2oo, allowing the gas to escape from the 



gravitational influence of the halo. At the largest scales 
after ~ 3 Gyrs, we obtain converged values for the wind 
velocities of v ~ 45 km s _1 and 135 km s _1 for our two 
haloes (10 9 h' 1 M© and 10 10 /i _1 M©), respectively. (Note 
however that owing to our simplified initial conditions, our 
simulations cannot be taken as a realistic picture for radii 
larger than a few R200 where the wind will encounter the 
anisotropic mass distribution of the cosmic web.) 

Comparing the CR-streaming simulations of the dif- 
ferent halo masses, we observe a quasi-spherical outflow in 
the dwarf halo with 10 9 h~ x Mq. In contrast, there is a bi- 
conical, hour-glass shaped wind in the 10 10 /i _1 M© halo 
within R200 that starts to become more spherical on larger 
scales. This different wind morphology can be traced back 
to the density distribution of the star-forming ISM which 
determines the initial conditions at the base of the outflow. 
While the shallow potential of the dwarf halo is not able 
to confine the ISM into a thin disc and instead has an ISM 
with a quasi-spherical morphology, the larger halo shows a 
well-confined galactic disc. Such a configuration has a verti- 
cal stratification of CR pressure which determines the initial 
direction of the CR streaming and focuses the forming wind 
into azimuthal symmetry. 

We now detail the physics of the launching mechanism 
of the wind. Initially, the CR pressure distribution peaks at 
the star forming regions in the galactic disc which accelerate 
them in SN remnant shocks. This causes a steep CR pres- 
sure gradient that implies the onset of CR streaming out 
of the galaxy. A fluid element above the disc which is in 
approximate hydrostatic equilibrium with the surrounding 
fluid elements receives the (instreaming) CR flux and be- 
comes overpressured. CR loss processes (wave heating as a 
result of CR streaming instabilities or particle-particle inter- 
actions) transform a part of the CR pressure to the thermal 
plasma. Since the thermal plasma has a harder equation of 
state, it gains more in terms of pressure than the CRs have 
lost by this transfer (see Section [2J . This increases the to- 
tal pressure of this fluid element furthermore. This excess 
pressure causes a weak shock that accelerates the gas in the 
stratified atmosphere. Momentum conservation ensures that 
this weak shock moves in the same (outwards) direction as 
the original CR-streaming flow. The rarefaction wave, mov- 
ing in the opposite direction, dilutes the gas in the tail and 
causes a buoyant rearrangement of the flow. We see that the 
streaming CRs are crucial in this picture for launching the 
outflow. 

A central question for winds is the evolution of their 
outflow velocity along the streamlines and whether they are 
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Figure 3. Maps of the baryonic gas density and the gas velocity field (indicated by arrows) for a halo of total mass 10 9 /i _1 Mq (top 
panels) and 10 10 /i _1 Mq (bottom panels) at t = 6/i -1 Gyr. The maps show an edge-on slice through the midplane of the galaxy that 
forms in the centre. The left panels correspond to the simulations with CR streaming and advection, while the right panels show the 
simulations with purely advective CR transport. While slow (thermal) gas motions dominate the velocity field in the advection-only run, 
the simulation employing CR streaming and advection shows a powerful outflow. 



energy or momentum driven. In the first case, it is known 
to be difficult to obtain the observed large mass loading 
factors. Additionally, if the (uncertain) clumping factor in 
the wind is larger than unity, the wind would radiate away 
more energy per unit time than in the smooth case which 
could potentially stall it. In Fig. 2] we show the radial ve- 
locity of the gas in a cone of opening angle 20°, aligned 
with the angular momentum axis of the disc. In the two 
lower-mass haloes, we see an accelerating wind over the en- 
tire radial range. This is due to the particular property of 
CR streaming-driven winds, which represents neither a clas- 
sic energy-driven nor a momentum-driven wind in a sense 
that they are not preloaded with those conserved quanti- 
ties. Instead, they are continuously re-loaded with energy 
and momentum during their ascent in the gravitational po- 
tential through the Alfven-wave heating term in the energy 
equation and the VP cr term in the momentum equation. 

While there was initially a wind launched in the 
lO n /i -1 M0 halo that propagates on scales r > 0.7i?2oo 
(at 2.5 h~ x Gyr), it apparently stalled for r < 0.7i?2oo and 



started to fall back as a fountain flow onto the galactic disc 
(see Fig. [4}. We observe the same behaviour for later times 
in the 10 10 /i _1 M© halo. The reason for this is the increasing 
stellar mass and the formation of a thin galactic disc which 
causes the wind to be launched deeper in the gravitational 
potential (implying a small r&) and to loose more binding 
energy during its ascent in the potentialo Another effect 
that can partly affect the sustainability of winds in larger 
haloes are the increased CR losses due to Alfven-wave heat- 
ing in the dense ISM which is characterized by steep CR 
pressure gradients. For our current description of the ISM 
with an effective equation of state, we may overestimate the 
CR losses as CRs are kept for a longer time in the dense 
phase of the ISM. In reality, they may be able to escape 



2 In both cases, the thermal pressure distribution is monotoni- 
cally decreasing for increasing r so that there is no pressure barrier 
building up in the halo due to CR Alfven-wave heating that could 
potentially stall the wind. 
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Figure 4. Radial velocity of the gas in a cone of opening an- 
gle 20°, aligned with the angular momentum axis of the disc 
for our three haloes of masses 10 9 /i -1 Mq, 10 10 /i -1 M@, and 
1O 1:l /i _1 M0 at 2.5/i _1 Gyr. In the two lower-mass haloes (as 
well as for r > 0.7i?200 in the 10 11 h~ 1 M@ halo), we see an 
accelerating wind over the entire radial range which is due to a 
continuous CR momentum and energy deposition during the as- 
cent of the wind in the gravitational potential. For r < 0.7i?200 
in the 10 11 h~ x Mq halo, the wind stalls and falls back onto the 
disc (shown by the negative velocity values) which is the charac- 
teristics of a fountain flow. 



more quickly into the warm and more dilute phase of the 
ISM that is characterized by a smoother CR pressure gra- 
dient so that less CR energy is lost at the initial phases 
of the wind launching. Without improving our modelling of 
the ISM or the CR transport within the ISM the resulting 
wind velocities could thus be underestimated. We will re- 
turn to this point in Section 15.21 and estimate its effect on 
the cumulative mass loss in a 10 9 h~ x Mq halo. 

To highlight the launching mechanism of the wind, we 
study the CR-to-thermal pressure ratio, X CT = P cv /P, in an 
edge-on slice through the galactic disc in Fig. \5\ Generally, 
we find an increasing CR pressure fraction, X cr , at larger 
radii. As the wind propagates, the composite gas of CRs and 
thermal plasma experiences adiabatic expansion so that the 
pressure drops less quickly than the thermal pressure due 
to the composite's softer equation of state. Thus, the CR- 
to-thermal pressure ratio rises to X cr > 1 at large radii. 
Apparently, this effect wins over CR energy losses due to 
CR Alfven-wave heating during the ascent of the wind in 
the halo potential. 

The X cr map of the 10 9 Ii~ 1 Mq halo shows a ho- 
mogeneous morphology for r < R200 with values around 
X cr c± 0.5. At larger radii, the morphology of the X cr map 
becomes patchier. It is interesting that CR streaming is 
unable to smooth this inhomogeneous CR pressure distri- 
bution. This is because of the large wind velocities which 
reach values of v ~ 50 km s _1 , whereas the sound speed 
(equal to the streaming speed in our model) is only around 
c s ~ (5 — 10) km s _1 or less. Thus, advection dominates the 
transport on these scales which is not expected to result in 
a smooth distribution of CR pressure and apparently does 
not significantly mix. The X CY map of the 10 10 h~ x Mq halo 
shows the hour-glass morphology and supports the picture 
of a bi-conical outflow driven by CR streaming. Clearly vis- 



ible are the torus-shaped regions in blue on scales below the 
virial radius which indicate a low relative CR pressure. This 
is due to a converging vortex flow towards the midplane and 
a necessary consequence of the bipolar outflow. Here, the 
CR pressure component is disfavoured in comparison to the 
thermal component upon adiabatic compression. 

5.2 Mass loss due to galactic winds 

After studying the launching mechanism of galactic winds 
through CR streaming, we will now quantify the associated 
mass loss and the mass loading of the wind. In Fig. [6] we 
show the fractional mass as a function of time for a refer- 
ence simulation without CR feedback, for a simulation with 
advective CR transport-only, and for two runs that addition- 
ally include CR streaming and differ only in the adopted CR 
acceleration efficiencies of Csn = 0.1 and 0.3. We define frac- 
tional mass as the baryonic mass in form of stars and gas 
contained within a certain radius at a given time, normal- 
ized to the total gas mass within this radius at the beginning 
of the simulation. 

The reference run without CR feedback exhibits no 
mass loss and so does the ad vect ion-only run, in accordance 
with the findings in Section 15.11 The fact that there is no 
wind in the former is connected to the implementation of 
stellar feedback that we use, which does not allow for out- 
flows without additional adjustments. In contrast, we de- 
tected a strong wind in our simulation with CR streaming 
which results in a considerable mass loss (see Fig. [6]). De- 
pending on the assumed CR acceleration efficiency, a frac- 
tion of (40 — 60)% of the original gas mass contained within 
the virial radius of the 10 9 h~ x Mq halo at the start of the 
simulation is expelled until t ~ 6/i _1 Gyr. The non-zero 
slope of the mass loss history indicates a moderate contin- 
uation of the mass loss if we increased the run time of our 
simulation. 

The mass loss history shows a strong dependence on 
halo mass with a fractional mass loss of only ~ 5% from 
the virial radius in the case of the 10 10 and 10 11 h~ x Mq 
haloes (see right panel of Fig. [6}. Thus, these outflows are 
weaker in comparison to that in the dwarf halo of M200 = 
10 9 /i _1 Mq, and do not result in a severe mass loss. Since the 
mass loss history has saturated in the two high-mass haloes, 
further mass loss is not expected for increasing simulation 
time. Apparently, the gravitational attraction is too strong 
for the CRs to launch a strong wind in these cases, but 
certainly drives powerful fountain flows. 

However, the mass loss history in Fig. [6] demonstrates 
that the outflow is not able to develop until t ~ 1.8 h _1 Gyr 
which is due to the ram pressure from inflowing gas that 
the wind needs to overcome first. Ram pressure can there- 
fore significantly delay wind formation and is probably more 
important in this respect than the shallow gravitational 
potential of this 10 9 h~ x Mq halo. This result is in agree- 
ment with previous stu dies of outflows ([Fuiita et al.ll2004l : 
iDubois & Tevssierll2008r ) and may be partly due to our sim- 
plified initial conditions. If most of the gas mass is accreted 
onto a galaxy along filaments that cover a small solid an- 
gle, the ram pressure in directions other than the filaments 
would be substantially reduced, possibly allowing for an ear- 
lier onset of a wind after the first burst of star formation. 

It is interesting to compare the mass loss rate in the 
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Figure 5. CR-to-thermal pressure ratio, X cr = P cr /P, in an edge-on slice through the galactic disc. Shown is the simulation with CR 
streaming and advection in haloes of 10 9 h~ 1 Mq (left) and 10 10 h~ 1 Mq (right) at time t = 6h —1 Gyr. The CR-to-thermal pressure 
ratio is less than 50% in the vicinity of the centre, because of loss processes that effectively transfer CR energy into the thermal reservoir, 
but becomes dominant at larger heights due to the softer adiabatic index of CRs. 
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Figure 6. Time-evolution of the baryonic mass (in form of gas and stars) enclosed within the virial radius, M(t), normalized to the 
initial mass at the beginning of the simulations. Left: we compare M(t) for four different scenarios (as indicated in the legends) in a halo 
of total mass 10 9 /i _1 M . Right: we show M(t) for three different halo masses, 10 9 h' 1 M , 10 10 h' 1 M , and 10 11 h~ x M , always 
adopting the CR streaming and advection scenario with £sn = 0.3. 
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Figure 7. Mass loss rates and SFRs as a function of time for the CR streaming and advection case, computed within the virial radius 
in haloes of total mass 10 9 /i _1 M (left) and 10 10 /i _1 M (right). The mass loss rate was calculated for two different acceleration 
efficiencies, £sn = 0.1 and 0.3. 



© 2008 RAS, MNRAS 000.ITU251 



Galactic winds driven by cosmic-ray streaming 11 



wind with the SFR. These are expected to be roughly pro- 
portional to each other, as simple theoretical predictions and 
observations suggest. Thus, in Fig. [T) we show the time evo- 
lution of the SFR and the mass loss rate within the virial 
radius for two CR streaming simulations that differ only 
in the acceleration efficiency of CRs. The influence of ram 
pressure is apparent, which delays a significant escape of gas 
until t ~ 1.8 /i -1 Gyrjj Once the CR-driven wind overcomes 
the ram pressure, the mass loss rate increases sharply. 

In case of the high acceleration-efficiency run for the 
dwarf haloes of 10 9 h~ x M©, the mass loss rate is about 5 
times the SFR and about 1.5 times the SFR for the lower 
acceleration efficiency (we always consider the "peak mass 
loading", i.e., we divide the respective maxima of the mass 
loss and star formation histories). For dwarfs, the mass load- 
ing of these winds is therefore considerable. In the end of our 
high acceleration-efficiency simulation, we find that the wind 
has expelled a gas mass fraction of ~ 60% while c± 15% of 
the initial gas mass ended up in stars (which was consid- 
erably smaller in comparison to the stellar mass fraction of 
~ 70% for our run without CR feedback, see Fig. [8}. In 
contrast, the wind in the 10 10 /i _1 M halo only has a small 
mass loading of about 0.13 for a high acceleration efficiency. 
While only ~ 5% of the initial gas mass was expelled by the 
wind in this halo, the stellar mass fraction dropped from 
~ 80% down to ~ 60% due to feedback by CR streaming 
as well as by advected CRs (Fig. [8}. Apart from that, in 
the low acceleration-efficiency CR streaming run there is no 
mass loss out of the virial radius. 

The energy transfer from the CRs to the thermal plasma 
following the excitation and damping of hydromagnetic 
waves within the dense ISM, that has a comparatively short 
cooling time, causes a large fraction of this energy to be lost 
to cooling radiation. If CRs are injected into the warm phase 
of the ISM from where they would drive a wind, then our 
simplified treatment of the multiphase structure of the ISM, 
which does not explicitly treat the warm and hot gas in su- 
perbubbles, underestimates the effect of CR streaming and 
hence the associated mass loss rates. To estimate the impor- 
tance of this effect, we switch off the conversion of CR to 
thermal pressure in our CR streaming implementation and 
find a fraction of 90% of the original gas mass contained 
within the virial radius of the 1O 9 /i _1 M halo is already 
expelled by t ~ 4/i _1 Gyr. For the reminder of this work, 
we always include CR-wave heating but bear in mind that 
the results may be too conservative as a result of too effi- 
cient CR energy losses to radiation in the dense, cool phase 
of the ISM. 



5.3 Suppression of star formation rates 

To study how CR feedback influences the star formation pro- 
cess, we plot the star formation history in Fig. [8] We vary the 



3 In both simulations, there is already a small amount of gas 
leaving the halo around t ~ 0.7 h~ 1 Gyr. This is not due to a wind 
but rather caused by the build-up of thermal and CR pressure in 
the centre of the galaxy owing to star formation, that reaches its 
maximum at the same time. This will slowly push a small amount 
of gas over the boundary of the halo. However, the negative mass 
loss rate immediately thereafter indicates that this gas flow does 
not propagate far, but returns soon. 



simulated physics and compare a model without CR feed- 
back, one with advective CR transport only, and one that ad- 
ditionally includes CR streaming. For each of our CR mod- 
els, we employ two CR acceleration efficiencies of Csn = 0.1 
and 0.3, respectively. Obviously, if feedback by CRs is taken 
into account, the SFR is suppressed in comparison to the ref- 
erence simulation without CRs. In case of the CR streaming 
model with a high acceleration efficiency, the SFR is sup- 
pressed by about 80% (40%) at the peak of the star forma- 
tion history for our halo with 10 9 /i _1 M (10 10 h' 1 M ). In 
our model that only accounts for advective CR transport, 
this suppression is even increased and amounts to more than 
90% (70%) for our halo with 10 9 h' 1 M (10 10 / i " 1 M ); in 
accordance with the findings of ljubelgas et al.l (|2008). As 
laid out there in detail, the suppression of the SFR by CR 
feedback is less effective in higher mass haloes which attain 
higher central gas densities due to their deeper potential 
wells. This is because adiabatically compressing a composite 
of CRs and thermal gas disfavours the CR pressure relative 
to the thermal pressure due to the softer equation of state of 
CRs. Hence, for the higher central densities present in larger 
haloes the relative importance of the CR pressure as well as 
their modulating effect on the SFR decreases. 

There are two reasons why the inclusion of CR stream- 
ing suppresses the SFR in comparison to the reference model 
without CRs. First, CR pressure that quickly builds up as 
stars form and eventually explode will hinder the gas from 
collapsing to densities much larger than the threshold of star 
formation. This is due to the inefficient cooling processes 
of CRs as opposed to the thermal gas that cools compara- 
tively quickly. Consequently, the SFR, that scales with the 
gas surface density as Ssfr oc Sg a t, will be reduced. As 
soon as CRs are removed by cooling processes, collapse can 
start again and the SFR increases. The overall effect is a 
suppression of star formation as well as an oscillatory be- 
haviour of the SFR, that can be appreciated in Fig. [8] If CR 
transport via streaming is included, these oscillations basi- 
cally vanish because freshly accelerated CRs continuously 
stream away from the star forming regions, allowing the gas 
to reach higher densities and consequently enhancing the 
SFR in comparison to the advect ion- only run. 

Second, the wind launched by CR streaming prevents 
a large fraction of the infalling gas from accreting onto the 
galactic disc. Hence, the amount of gas accumulated in the 
galactic disc is reduced. This lowers the available amount 
of gas that may be turned into stars as well as the surface 
gravity. The resulting smaller central gas densities imply a 
lower SFR. 

The runs adopting a low CR acceleration efficiency of 
Csn = 0.1 for the CRs show less suppression of the SFR in 
comparison to their corresponding high-efficiency counter- 
parts. This is due to the lower contribution of CR pressure 
to the total pressure in the star forming regions, weakening 
the influence of CRs on the process of star formation, as 
outlined above. Additionally, the mass loading of the winds 
is also reduced in this case. Nevertheless, the suppression is 
still significant and amounts to ~ 40% (~ 80%) at the max- 
imum of the star formation history in the CR streaming 
model (advection-only model) for the 10 9 h~ x M halo. 

The bottom panels of Fig. [8] show the time-evolution of 
the total amount of mass lost from the halo's virial radius 
for our CR streaming simulation with £sn = 0.3 as well as 
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Figure 8. Top panels: SFR as a function of time in a halo of total mass 10 9 /i -1 Mq (left) and 10 10 /i _1 Mq (right). We compare five 
different simulations: a reference simulation without CRs (black line) and four runs with CR feedback, differing only in the acceleration 
efficiency and the transport scheme, i.e., advective transport-only (orange line: £sn = 0.1, green line: £sn = 0-3) and CR streaming and 
advection (blue line: Csn = 0.1, red line: Csn = 0.3), respectively. CR feedback leads to a strong suppression of the SFR in comparison 
to the reference simulation. Bottom panels: Time-evolution of the total mass lost in the wind out of the virial radius R200 (dashed red 
line), and of the total mass in stars (solid red line) for the CR streaming run. These are compared to the total mass in stars in the 
advection-only run (green line) and in the reference simulation without CRs (black line). Both runs with CRs adopt £sn — 0.3. In the 
10 9 h~ x Mq halo, the amount of mass removed from the halo in the CR streaming run is considerably larger than the stellar mass that 
forms in this simulation, in agreement with the large mass loading of the outflow. Moreover, the total mass in stars is significantly less 
than in the reference simulation. 



the total amount of mass formed in stars at a given time in 
this simulation. We contrast this to the cumulative stellar 
mass for a model without CRs and our model with advective 
CR transport-only. This shows again that the mass loading 
of the outflow in the CR streaming case is large, with the 
total mass loss being about 4 times as large as the mass in 
stars at the end of the simulation. Furthermore, the final 
stellar mass at the end of the CR feedback runs is reduced 
by ~ (80 — 90)% relative to the final stellar mass in the ref- 
erence simulation (considering the 1O 9 /i _1 M0 halo). This 
is in good agreement with the suppression factor of the SFR 
at the peak of star formation history. As expected, the re- 
duction is strongest in the model that exhibits only advec- 
tive CR transport which results in the least amount of stars 
formed at the end of the simulation. 

In contrast, the baryonic mass lost in the galactic wind 
of the larger halo of 10 10 h' 1 M is only - 8% of the total 
mass in stars at t ~ 6/i _1 Gyr, reflecting the small mass 
loading of the wind (see Section T5.2[) . In this larger halo, the 
final stellar mass at the end of both CR feedback runs is 
only reduced by ~ 20% relative to the final stellar mass in 
the reference simulation. The large reduction of the SFR at 



the peak of ~ 70% in the advection-only run is therefore not 
seen in the final stellar mass. This is because of the different 
star formation history which shows more continuous rather 
than a bursty behaviour of the CR streaming models for this 
halo. 



5.4 How CR streaming heats the halo gas 

In this section, we assess the role of the wave heating that 
is intimately connected to the physics of CR streaming. To 
this end, we only concentrate on our CR streaming model 
that successfully launches winds. In Fig. [9] we plot the gas 
temperature in an edge-on slice at two different epochs, t — 
2.8 h' 1 Gyr and t = 3.2 h' 1 Gyr, about 1 h' 1 Gyr after the 
outflow started. At t — 3.2 /i -1 Gyr we can see collimated 
flows through chimneys of hot gas, extending up to \z\ ~ 
20/i -1 kpc above and below the disc. In these chimneys, the 
gas temperatures reach values greater than T ~ 10 5 K. In 
the centre of the maps, we can see that the disc is also heated 
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Figure 9. Heating mechanism of the halo gas in a halo of total mass 10 10 h~ x Mq for the CR streaming model. Top panels: Time evolution 
of the gas temperature and the gas velocity field in an edge-on slice through the galactic disc. Results are shown at two different epochs 
in the simulation, at t = 2.8 h —1 Gyr (left panels) and t = 3.2 /i —1 Gyr (right panels). The maximum velocities encountered in these runs 
are 26 km s _1 (left panels) and 70 km s _1 (right panels). Collimated cavities of hot gas form below and above the disc and propagate 
further up in the course of time. Middle panels: Corresponding plots for the time-evolution of the ratio of the wave heating rate due 
to the damping of self-excited waves to the gas cooling rate, A waves /A coo i, demonstrating that wave heating is responsible for heating 
the cavities. Bottom panels: Gas density in a cutting plane parallel to the disc at a height of z = 4.5 /i -1 kpc (to assess the density 
distribution in the hot cavities). At t = 3.2 /i —1 Gyr, an under-dense channel forms in the centre where the cooling rate is consequently 
lowered so that wave heating can dominate over the reduced cooling rate. 
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Figure 10. Temperature distribution at the time of maximum CR Alfven-wave heating in an edge-on slice through the galactic disc. 
We compare three different haloes of mass 10 9 h~ 1 Mq, 10 10 h~ 1 Mq, and 10 11 h~ 1 Mq (left to right) in our CR streaming model with 
an acceleration efficiency of Csn = 0-3. Note that the resulting halo temperatures roughly scale as kT oc v^- d ~ v esc- The temperature 
structure resembles that of the wind, which implies that with increasing halo size, the morphology of the hot patches becomes more 
conical. The broadening of the hot regions in our 10 11 h~ 1 Mq halo is associated with the inability of CR streaming to drive a sustained 
wind that escapes from such a halo. Hence, the kinetic energy of the fountain flow drives turbulence which dissipates energy and thereby 
heats larger regions of the halo gas. 



up to T ~ 10 5 K owing to supernova feedback |j These hot 
regions are surrounded by a torus-shaped region, filled with 
comparatively cool gas at T ~ 10 4 K. 

The only conceivable process which could provide a 
source of heating several kiloparsecs above the disc, where 
the gas density is low, must be the damping of CR-excited 
waves. In Sec. [2] we have seen that the hydromagnetic waves 
are created by streaming CRs in the background plasma. 
Those waves are damped on short time-scales, thereby trans- 
ferring energy to the plasma. This heating process is in oper- 
ation as long as CRs are streaming through the rest frame of 
the gas, i.e., as long as there is a gradient in the CR pressure. 
Thus, unlike Coulomb or hadronic losses, which scale with 
gas density, wave heating is also important in low-density 
regions. In order to demonstrate that this process is respon- 
sible for the hot gas chimneys, we show an edge-on slice 
through the galaxy with the ratio of the radiative cooling 
rate of the gas to the wave-heating rate of equation (|A31[) 
due to the damping of self-excited waves, A W aves/A CO oi (mid- 
dle panels of Fig. [9}. 

At t = 2.8/i _1 Gyr, wave heating is nowhere able to 
overcome the radiative cooling of the gas. This is in partic- 
ular true for the galactic disc, where the density is so large 
that the cooling dominates the heating by more than two or- 
ders of magnitudes. At t = 3.2 /i _1 Gyr, however, we can see 
the structure of the hot gas chimneys again, traced by the 
evidence that the wave heating dominates the cooling there 
by about a factor of ten, showing that heating via wave 
damping is indeed the process that creates the hot cavities 
in our simulations. 



4 W e are using a sub-resolution m odel of the ISM and star forma- 
tion (|Springel fc Hernquistll2003r ). It assumes that star formation 
establishes a self-regulated regime of the ISM, which arises due to 
the interplay of cold, dense clouds that continuously form stars 
and supernova feedback that evaporates these clouds, creating a 
hot ambient medium. The sub-resolution model describes the two 
phases by an effective equation of state that is characterized by 
a mean mass weighted temperature which we plot in Fig. [9] 



To better illustrate the reasons for the dominance of 
wave heating over gas cooling above and below the disc, 
we take a slice parallel to the disc that cuts the outflowing 
gas stream at a height of z = 4.5 /i -1 kpc above the disc, 
which is near the base of the upper chimney. We plot the 
gas density in this slice in the bottom panels of Fig. [9] At 
t = 2.8 /i _1 Gyr, when the cavities have not yet formed, we 
see that the outflow covers a circular area in the plane and 
features a denser core in the centre, with decreasing density 
towards the outskirts. This indicates that the outflow is very 
collimated. Later, at t — 3.2 /i -1 Gyr, an under-dense hole 
has appeared in the very centre of the core. Thus, the gas 
cooling rate will drop there, so that the CR-wave heating 
can heat up the gas. 

What is the reason for this under-dense channel in 
the centre of the outflowing material? Interestingly, the 
hot chimneys first occur about l/i -1 Gyr after the outflow 
started. The likely reason for this is that the CR-driven wind 
becomes more and more collimated as time progresses, ow- 
ing to the disc that forms in the centre of the simulation box. 
When the wind first occurs, the disc formation is not yet fin- 
ished, so that it still has an approximately spherical shape. 
At these early times, CRs can stream at a large angle with 
respect to the z-axis without encountering too much disc 
material, resulting in an outflow with a wide opening angle. 
Later, when the disc has formed, the inertia of the dense and 
cool star forming gas prohibits a long pathway of the CRs 
through the disc and forces an outflow with a smaller open- 
ing angle. When the outflow is collimated enough, it will 
then quickly dig a diluted channel in the old eject a above 
the disc with a correspondingly lower cooling rate. 

This argument is also in agreement with the fact that 
we do not observe chimneys of hot gas for our lower mass 
halo. The weaker gravity of those dwarf haloes does not sup- 
port a thin disc forming at the centre, but a rather spherical 
density distribution which is unable to collimate the wind. 
This is highlighted in Fig.[10]which shows temperature maps 
for haloes with 10 9 /i" 1 M , 10 10 /i" 1 M , and 10 11 /i" 1 M . 
These show a clear dependence of the wave heating on the 
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halo mass. The resulting halo temperatures roughly scale 
as kT oc ^wind ~ ^esc with the largest halo in our simula- 
tions reaching temperatures in excess of 10 6 K so that the 
outflow is expected to emit thermal bremsstrahlung emis- 
sion. Besides, the larger SFR obtained for the higher mass 
halo and correspondingly the higher CR energy densities 
could contribute to the preferred creation of the chimneys 
in the higher mass halo. We also see in our largest halo 
that the heated regions are less collimated in comparison 
to the intermediate-mass haloes: the wind starts deeper in 
the gravitational potential of this halo (see Section [3]) and 
looses a good fraction of its kinetic energy in climbing up the 
greater potential difference so that it is unable to counteract 
the ram pressure of the infalling gas (that reaches also bigger 
infall velocities in assuming the larger gravitational binding 
energy of this halo). However, the interaction of the outflow 
with the infalling gas results in violent turbulence that also 
dissipates by cascading down in length scale, thereby ad- 
ditionally heating the halo gas. Potentially, the additional 
momentum deposition from radiation pressure may help in 
circumventing the stalling of the CR-driven wind seen for 
this halo. 

Another process that certainly enhances the collima- 
tion of the wind into a narrow channel is directly connected 
to the physics of CR streaming. Recall that the magnitude 
of wave heating due to CR-excited hydromagnetic waves is 
proportional to the streaming speed which is equal to the 
sound speed in our model. Thus, CRs will stream fastest 
and deposit most of their energy in the high-temperature 
gas inside the channel, thereby further increasing the tem- 
perature there and more importantly the outflow velocity 
because of enhanced thermal pressure. This is illustrated 
by the increasing magnitudes of the velocity vectors in the 
outflow (see top panels of Fig. [9J. We find maximum veloc- 
ities of rsj 70 km s _1 inside the hot channel, providing the 
gas with sufficient momentum to overcome the gravitational 
attraction of the halo. 

Moreover, one can see vortex like structures in the ve- 
locity field at the sides of the chimneys, where gas is stripped 
off from the flow and returns to the disc. The vortical mo- 
tions are presumably created by shear along the sides of the 
outflow and they converge towards the base of the cavities, 
where gas is then compressed to high densities, facilitating 
efficient radiative cooling (similar effects a lso occur in the 
case of hot, rising supernova bubbles, see iRobinson et al.l 
1 20041 ). This converging gas motion therefore explains the 
cool torus structure surrounding the disc. As time pro- 
gresses, the vortical motions will rip the warm chimneys 
more and more apart so that a warm shell of gas is created 
around torus and the galaxy. 

5.5 Maps of Ha emission 

In the previous section we have seen that the wave heat- 
ing owing to CR streaming results in characteristic, high- 
temperature structures. In these structures, the hydrogen 
gas is highly ionized and therefore allows for optical line 
emission of recombining electrons. Important in this context 
is the Ha line (the dominant line of the Balmer series cor- 
responding to the 3—^2 transition in the hydrogen atom), 
because hydrogen is the most abundant element in the Uni- 
verse and thus ubiquitous in astronomical spectra. Ha-line 
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Figure 11. Ha emission (edge-on slice) in our CR streaming sim- 
ulation at t = 4.5 h~ 1 Gyr for a total halo mass of 10 10 /i -1 Mq. 
We see diffuse Ha emission in the halo (R200 ~ 35 /i -1 kpc), 
which is explained by ionized gas, entrained in the wind. 



emission from hot gas has been detected in many galaxies 
that ex hibit a wind, e.g., in the well-known starburst galaxy 
M82 (JBland fc Tullvlll988h . Thus, we calculate a map the 
expected Ha emission in our CR streaming and advection 
simulation to see how well they compare to these observa- 
tions. 

The Ha emissivity, e (Ha), is given by (|Spitz 



~3r1ll978h 



e (Ha) = hu a a eS (T) n e n p 



(17) 



where hv a = 1.88 eV is the energy of a Ha photon, a e s(T) is 
the effective recombination coefficient and n e n p denotes the 
product of the (free) electron number density and the pro- 
ton number density, respectively. The product a e ff(T)n e n p 
gives the number of Ha photons emitted per second and 
per cubic centimetre. The temperature dependence of the 
effective recombination coefficient can be fitted, 



a eff (T) = 10 



_i3 2.708 T 4 -°- 648 



3 -1 

cm s 



(18) 



1 + 1.315 T 4 0523 

with T 4 = T/10 4 K (|Peauignot et al.lll99lL assuming case 
B). Using equation (|18[) we can then integrate the emissivity, 
equation (|17p . along the line- of- sight. 

The result is shown in Fig. 111! where we show the 
Ha emissivity in an edge-on view of our simulation at 
t — 4.5 h~ x Gyr, when the wind is already well-developed. 
The disc dominates the surface brightness map due to the 
supernova feedback which heats the gas to sufficiently high 
temperatures so that the gas becomes ionized. As a result, 
the product n e n p is large and consequently also the recombi- 
nation rate. Interestingly, there is a halo of Ha emission with 
a characteristic surface brightness of ~ 10 -13 erg cm -2 s _1 
that extends to scales of ~ 20h~ 1 kpc. Such a diffuse emis- 
sion might be detect able, given reported detections at com- 
parably low levels (| Young et al.l 1 19961 ). This emission is 
caused by the outflow that transports ionized gas into the 
halo and beyond. Especially inside the cavities and in their 
surroundings, where the streaming CRs heat the gas via 
wave excitation, the gas has a high degree of ionization, 
ranging from ~ 70% to up to 100% in some places. Unfortu- 
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nately, the overall gas density in the chimneys is very low, 
so that their clear structure that we saw in the tempera- 
ture maps is not visible in Ha. There is in addition even 
lower surface brightness diffuse emission with an elliptical 
morphology as a result of the galactic wind, extending up 
to \z\ -80/i _1 kpc. 

5.6 Limitations of the model 

We need to caution that some of our conclusions depend on 
certain assumptions used in this work and may not be fully 
generic to the physics of winds driven by CR streaming. (1) 
The sub-resolution model for star formation that we adopt 
here is arguably not well-suited for studying outflows. In 
particular, the scheme is designed to produce stars in a qui- 
escent mode, i.e., star formation in general only takes place 
at an average rate. In reality, how ever, winds are typi cally 
observed for starburst galaxies (cf. IVeilleux et al.ll2005l . and 
references therein), which exhibit high rates of star forma- 
tion and therefore inject a larger amount of thermal and CR 
pressure. Apart from that, the star formation scheme does 
not allow for an explicit treatment of supernova remnants 
and of the hot, over-pressurized gas bubbles they contain. 
These bubbles are expected to rise buoyantly in the disc 
potential of their host galaxy which should open channels 
along which the wind can escape much easier. (2) The use 
of such a sub-resolution model with an effective equation 
of state of the ISM may overestimate the energy transfer 
from the CRs to the thermal plasma following the excita- 
tion and damping of hydromagnetic waves within the dense 
ISM, that has a comparatively short cooling time, causing a 
large fraction of this energy to be lost to cooling radiation (as 
shown in Section 15. 2ft . This overestimate of CR streaming 
losses may be particularly important in larger haloes with 
M 2 qq > 10 10 M where larger ISM overdensities are reached 
due to the larger gravity. When explicitly modelling these 
two phases, it has been shown that CRs do not effectively 
couple to gas within cool clouds since they do n ot exert 
forces inside of cool clouds ([Everett fc ZweibefeOllh . Hence, 
if CRs were injected into the warm phase of the ISM (for an 
improved ISM modelling in the simulations), CR streaming 
could more efficiently drive winds, possibly producing larger 
mass loading factors for haloes with M200 > 1O 1O M0. (3) 
To isolate the physics of CR streaming, we have here used 
an extremely simplified model for a forming and evolving 
galaxy. Our initial conditions assume spherical symmetry 
which causes the gas accretion and the associated ram pres- 
sure to be initially spherically symmetric. However, a galaxy 
that forms in a cosmological environment accretes the gas in 
part along filaments that have a much smaller angular cov- 
ering factor, hence implying a ram pressure that varies with 
angle and could potentially modify the resulting wind mor- 
phology. This effect is particularly important for collimated 
winds in larger galaxy haloes with M200 > 10 10 M©. Those 
should feel a substantially reduced amount of ram pres- 
sure in the wind direction since filaments are generally not 
aligned with the disc's angular momentum axis. (4) In our 
approach, we use a simplified treatment of the CR stream- 
ing speed and equate it to the local sound speed. However, 
a realistic value of v s t should depend on the pre-existing 
wave level, thermal and CR energy density, magnetic topol- 
ogy, and the damping mechanism of the Alfven waves, to 



name a few. (5) We do not follow MHD; in particular, we 
assume the existence of locally open field lines along which 
CRs can escape into the halo. While only 10% of the galac- 
tic SN remnants will create flux tubes that reach heights 
> 10 kpc above the disc, SN remnant lobes are expected to 
overlap, suggesting that at every location there is su fficient 
supply of open field lines (|Breitschwerdt et al.lll993T ). This 
lets this assumption appear plausible (see also Appendix lA]) . 
(6) We neglect diffusive shock acceleration of CRs in winds 
from high-mass Wolf-Rayet stars. The recent detection of 
TeV-gamma rays from the young stellar cluster Westerlund 
2 is compatible with the hypothesis of pion-decay emis- 
sion of relativistic nuclei interacting with a close-by molecu- 
lar cloud complex (|HESS Collaboration et al.ll201lh . If con- 
firmed, that would imply that CR feedback could already 
be an important mechanism in dispersing molecular clouds 
when the most massive stars turn on, long before the first 
star became a supernova, making it possibly a competing 
mechanism to radiation pressure in launching outbursts. 



6 DISCUSSION AND CONCLUSIONS 

6.1 Properties of winds driven by CR streaming 

Using numerical simulations and analytical arguments, we 
show in this paper that CR feedback is able to power large- 
scale galactic outflows. These are not only able to escape the 
galactic disc, but also its host dark matter halo. In particu- 
lar, we demonstrate that CR streaming relative to the rest 
frame of the gas is necessary to drive winds through CRs. 
Solely considering advective CR transport, i.e., employing 
the flux-freezing approximation would not result in galactic 
outflows. In this advective-only case, the CR pressure does 
volume work on the gaseous discs (in the low-mass haloes), 
causing them to become more dilute which increases the 
thermal cooling time. As a result, the SFR declines strongly. 
After some time, the CR pressure is dissipated such that the 
gas can collapse again, triggering another "cycle" of star for- 
mation and causing an oscillatory mode of star formation 
with substantial suppression of the stellar m ass formed by a 
facto r of 8 for a halo of mass M200 < 10 9 M ([jubelgas et al.l 
2008). In this case, most of the CR energy is eventually 
transferred to the thermal gas, which radiates it away. In 
contrast, allowing for CR streaming in the rest frame of 
the gas down their pressure gradient enables CRs to quickly 
escape their acceleration sites in the dense, star forming re- 
gions, and to move into the CR-dilute, low density gas above 
the disc (or in between the cold phase of the ISM). There, 
the CR cooling rates (due to Coulomb or hadronic interac- 
tions) are significantly reduced since all these processes scale 
with the gas density. Thus, a big fraction of their initial en- 
ergy can be used to power the wind. Moreover, launching a 
wind from a low- density gas is considerably easier than from 
the dense phase, because the smaller column density of the 
material that needs to be removed. Star formation is still 
suppressed in comparison to the case without CR feedback, 
but less in comparison to the ad vect ion-only case. Remark- 
ably, the star formation rate is much smoother and does not 
show the oscillatory behaviour discussed above. 

We find that CR streaming powers winds most effec- 
tively in small dwarf haloes where the wind speed exceeds 
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the escape velocity for host halo masses M200 < 10 11 M©. 
Below this mass threshold, we find that the wind velocity 
increases as a function of host halo mass. We can trace this 
behaviour back to the potential difference the wind has to 
climb until it can escape, and which is effectively set by 
the scale height of the disc. On dwarf scales, the disc scale 
height is increasing with circular velocity so that the sys- 
tem is approximately scale invariant. In larger galaxies with 
M200 > 10 11 M©, the disc height is set by the effective equa- 
tion of state of the ISM and roughly constant with halo 
mass. This breaks scale invariance, causes an increasing po- 
tential difference for larger haloes, and renders CR stream- 
ing not to be powerful enough to be solely responsible for 
galactic winds. While CR-driven winds are approximately 
spherical in dwarf haloes (M200 — 10 9 M©), the denser discs 
in larger haloes are able to shield the solid angles subtended 
by the disc and hence focus the wind into bi-conical outflows 
around the angular momentum axis of the disc. 

In a halo of mass M200 = 10 9 /i _1 M© , the relative bary- 
onic mass loss amounts to ~ 60% (~ 40%) of the initial gas 
mass contained inside the halo's virial radius, assuming a 
CR acceleration efficiency of (sn = 0.3 (Csn = 0.1). Simi- 
larly, the mass loading factor of the wind increases towards 
low- mass haloes, reaching values as great as 5 for a halo of 
M200 — 10 9 h~ 1 M(? ) . The reason for these increased mass 
losses towards smaller halo masses are the shallower gravi- 
tational potential wells of these systems. However, despite 
the inability of CR streaming to drive winds in large galaxies 
that escape the halo, we nevertheless find powerful fountain 
flows. While we expect a smaller impact of CR streaming to- 
wards larger halo masses M200 > 1O 1O M , we caution that 
our simplified simulation setting may artificially weaken the 
CR streaming efficiency of driving winds in large galaxies. 
First, our employed sub-resolution model with an effective 
equation of state of the ISM may overestimate the CR en- 
ergy losses in the ISM due to Alfven wave excitation since 
larger ISM overdensities are reached in larger galaxy haloes 
due to their greater potential depths. Second, a cosmologi- 
cally forming galaxy accretes the gas in part along filaments 
that have a much smaller angular covering factor. Since fila- 
ments are generally not aligned with the disc's angular mo- 
mentum axis, this reduces the accretion ram pressure that 
collimated winds (aligned with the angular momentum axis) 
have to exert pdV-work on. Future work that realistically 
models the ISM in a galaxy forming in a cosmological set- 
ting is required to quantify these effects. 

Another unique property of CR-driven winds is the ef- 
fective heating of the outflow by dissipating Alfven waves 
that are excited by streaming CRs. Because the wind speeds 
increase as ^ w ind oc ^ eS c, we also see an increased rate of 
heating with increasing halo mass that results in tempera- 
tures of the halo gas of kT oc t^ind- We predict outflows in 
haloes with masses M200 > 10 10 M to reach temperatures 
in excess of 3 x 10 5 K so that the outflow is expected to emit 
thermal bremsstrahlung emission. 

We would like to emphasize that neither the classic 
energy-driven nor the classic momentum-driven wind fully 
captures the physics of CR streaming-driven winds. The 
classical picture assumes a pre-loading of the wind with 
energy or momentum. Instead, the physical mechanism of 
CR-driven winds implies wave excitation by the CR stream- 
ing instability that transfers both, energy and momentum to 
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Figure 12. Contours of v w i n d f° r c °ld (T ~ 10 4 K) wind models 
driven by a combination of r adiation and ram pressure (equation 
(12) of S harma fc Nath[2012h . Different symbols indicate observa- 



tions of winds with (par t ly) considerable mass loadings (IMartin 
I l998| : [Genzel et al Il200ll: ISchwartz fe Martini 120041 : iRupke et al" 
l2005l : lMartinll2005l : IWeiner et al.ll2009h . In these models, the re- 
gion to the lower right of these contour lines allows for powerful 
winds driven by ram and radiation pressure, ejecting the gas from 
the halo. In the region to the upper left, these two processes are 
not able to drive galactic winds, leaving the observational data in 
this region unexplained. The range of SFRs encountered in our 
CR streaming simulations is shown with red lines and indicates 
the importance of this process in driving winds, especially for 
galaxies showing smaller values of v c and SFRs. 



the gas. This manifests itself through the continuous Alfven- 
wave heating term in the energy equation and the VP cr 
term in the momentum equation. CR-driven winds are a 
two-stream phenomenon composed of the gas flow (carrying 
momentum) and the fast CR stream (exerting pdV-work on 
the gas) which continuously re-powers the wind with energy 
and momentum during its ascent in the gravitational poten- 
tial. 



6.2 Comparing different physical mechanisms for 
galactic winds 

We now scrutinize different physical mechanisms responsi- 
ble for driving galactic winds and put CR streaming - the 
mechanism studied here - into the context of other plausible 
processes, namely winds driven by radiation or ram pressure. 
Generally, it is believed that galaxies have winds when their 
star formation per unit area exceeds a threshold value of 
0.1 M© yr _1 kpc -2 . In our simulation of CR-driven winds, 
the SFR densities are lower than this threshold. Therefore, 
this process can address winds with low SFRs. 

In Fig. 1121 we show contours of ^ w ind in the parameter 
space spanned by v c and SFR for cold (T ~ 10 4 K) wind 
models driven by a combination of ra m pressure and radi- 
ation pressure ([Sharma &; Nathll2012h . According to these 
models, the region of parameter space to the lower right of 
these contour lines allows for powerful winds driven by ram 
and radiation pressure that are able to eject the gas from the 
halo, i.e., it is easier to drive winds for a smaller gravitational 
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potential well (smaller v c ) and a greater total power input 
to the system (larger SFR). In contrast, in the region to the 
upper left, these two processes are not able to drive galactic 
winds and demonstrate the failure of ram and radiation pres- 
sure driven wind models to explain the wind observations i n 
starburst dwarfs (IMartinl [l998l : ISchwartz fe Martini l2004h . 
Our simulation results for winds driven by CR streaming 
are shown with red lines (indicating the range of SFRs) and 
highlight the importance of CR streaming in driving winds, 
especially for galaxies showing small values of v c and SFRs. 

Good examples are the winds observed in the dwarf 
starburst galaxies NGC 4214 (SFR = 0.17 M^yr" 1 and 
ticmax = 40kms~ 1 : IWalter et al.ll2Q0ll : IMartinl 120051 ) and 
NGC 156 9 (SFR = (0.16 - 0.4) M^yr" 1 ari d i; c , max = 
35 km s~ 1 ; iMartin et al]l2002l ; IStil fe Israe1l2002h . The latter 
shows direct evidence for a metal-enriched wind from a dwarf 
starburst galaxy with a mass loading factor of 9 that carries 
nearly all the metals ejected by the starburst. Both galax- 
ies populate the no-wind region considering only ram and 
radiation pressure. These dwarf starbursts are close to our 
10 10 /i -1 M© halo in the parameter space shown in Fig. 1121 
The CR-driven bi-conical outflow in this galaxy resembles 
the observed winds in NGC 4214 and NGC 1569, consistent 
with a similar physical origin. 

While wind velocities and masses can be reliably de- 
termined for neutral phases (as it is the case for the two 
examples discus sed), and fo r Ha-emitting warm ionized 
gas (as done in IMartinl 1 19981 ). the velocity of the hot gas 
(10 6 — 10 7 K) is more difficult to determine. However, it 
is believed to be higher than that of the neutral gas, with 
velocities of order its own thermal sound speed of several 
hundred to a thousand kilometres per second. How does a 
neutral cloud (or some of its material) get entrained into the 
wind, in particular for a CR- streaming driven wind? Such 
a wind is necessarily magnetized with contributions from 
open field lines (extending from the disc into the halo) , flux 
frozen disc fields, as well as self- amplified field by means of 
CR streaming. This magnetized super-Alfvenic and super- 
sonic wind (with respect to the gas in the cloud that exhibits 
comparably low temperatures and sound speeds) impacts 
the neutral cloud and forms a contact discontinuity. Inde- 
pendent of the degree of magnetization, the cloud sweeps 
up enough wind magnetic field at the interfa ce to the wind 
to bu i ld up a dynami c ally import ant sheath (iGregori et al 



2Q0d : iLvutikovl 120061: lAsai et all 120071 : iDursi fe Pfrommei 
20081 : IPfrommer fe Dursill20ld ) 



The field strength of this magnetic draping layer is 
set by a competition between 'ploughing up' and slipping 
around of field lines, maintaining a magnetic energy density 
in steady state that is of order the wind ram pressure seen 
by the cloud. This is an inherently three-dimensional prob- 
lem where the third dimension is essential in allowing for 
field lines to slip around the obstacle. This strongly magne- 
tized layer modifies the dynamics of the cloud, pote ntially 
supp ressing hydrodynamic instabilities and mixing (|Dursil 
120071 ). and changing the geometry of stripped material. Most 
importantly, magnetic draping accelerates the cloud into the 
direction of the wind through magnetic tension of draped 
field lines that are anchored in the hot phase of the wind, 
streaming already ahead of the cloud (|Dursi fe Pfrommerl 
1 20081 ). Note that these authors study the opposite problem 
of a cloud moving through a magnetized wind (which can be 



transformed by means of a Galilean transformation to our 
problem at hand) and find that the deceleration through 
magnetic tension is always more important, by a factor of 
~ 3.7, than for the case of highly turbulent (Reynolds num- 
ber Re ~ 1000) hydrodynamic drag. The terminal velocity 
that the cloud can reach should be a sizable fraction of the 
wind velocity, and may depend on magnetic geometry and 
the magnetization of the wind. 

This emerging picture suggests that different physical 
processes dominate the launching and powering of the wind, 
with a relative contribution that depends on the SFR and v c 
(JHopkins et al.ll201ll ; Istiarma fe Nathll2012h . For increasing 
values of SFR and v c , we seem to have first CR stream- 
ing, then ram pressure, and finally radiation pressure as 
the dominant agent in powering winds. We conclude that 
(1) CR-driven winds are most effective in dwarfs, ejecting 
a large fraction (~ 60%) of metal-enriched baryons into the 
IGM, and (2) CR Alfven-wave heating is especially efficient 
in heating a fountain flow in high-mass galaxies. While we 
expect that our results for this physical scenario are ro- 
bust irrespective of redshift or environment, we caution that 
much more work is needed to firmly confirm this picture, es- 
pecially given the limitations in modelling these processes 
at the present time (see Section |5J 



6.3 Cosmological and observational implications 
of CR-driven winds 

In hierarchical structure formation, dwarf galaxies represent 
the building blocks of larger galaxies. CR-driven winds ap- 
pear to be very efficient in removing baryons from small 
haloes, hence this process should be the dominant feedback 
mechanism during the formation of these dwarfs, i.e., shortly 
before and after reionization. We expect the successive merg- 
ers of dwarfs to be less gas rich than without this feedback 
mechanism and to form less stars at high redshift. As galaxy 
systems grow larger, the ejected gas may be accreted later, 
but, due to the modified thermal history and angular mo- 
mentum, may end up at different places within the larger 
galaxies. Hence, it is not inconceivable that CR-driven winds 
in small-scale haloes have an indirect effect up to scales of 
L* at the knee of the luminosity function; but this needs to 
be carefully studied in future work. 

Additionally, we expect our CR-driven winds to influ- 
ence the early history of metal enrichment during the epoch 
of dwarf formation because of their considerable mass load- 
ing factors. This important role of low-mass systems in the 
metal enrichment of the IGM is in line with expectations 
from other authors (JNath & Trenthaml Il997l ; iMadau et"aH 
l200ll : ISamui et al.ll201(£ ~~ 

What about specific observational predictions of CR- 
driven winds? The bi-conical structure of the hot gas cavi- 
ties that formed in the 10 10 h~ x Mq halo has some similar- 
ities with H oj-structures, e.g., in the well-known starburst 
galaxy M82 (JBland fc Tullvl Il988l) or with the dwarf star- 
burst galaxy NGC 1569 (JMartin et al.ll2002h . However, the 
bi-conical structure observed in the temperature field has 
no direct resemblance in our simulated Ha emission maps, 
due to the low density of the gas in these cavities. This 
is the same problem that commonly makes observations of 
winds difficult, which face the difficulty of detecting the 
weak emission from a wind that is in general hot and di- 
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lute (see IVeilleux et al.l 120051 . and references therein) . Our 
simulations are not capable of describing the full spectrum 
of fluid phenomena that are likely to be associated with 
these hot cavities. Fluid i nstabilities of Kelvin - Helmholtz 
and Ray leigh- Taylor type ([Robinson et al.l 12004 ) as well as 
thermal instabilities should be important because of the 
density contrast around the chimneys and the shear along 
their sides. In addition, there will be entrainment of denser 
clouds from the interstellar medium when the hot wind 
pushes on them ([Mac Low fe McCravlll988h . All these phe- 
nomena will influence the morphology of the cavities and 
could enhance the density locally, so that their signature 
might indeed be visible in Ha. However, we lack the neces- 
sary resolution to see these effects in our simulations and the 
SPH method may also significa ntly contribute to this limita- 
tion, as recent studies found (Siiack i et al.ll201ll ). Moreover, 
magnetic fields that are confined within the hot chimney 
gas and rise up along with it are of particular importance. 
They will support the chimneys against the hydrodynam- 
ical instabilities and shear which try to tear them apart 
(iRobinson et al I l2QQ4l: iDursil 1200?! ; iRuszkowski et alJ l2007t 
iDursi & Pfrommerll 20081 ). Nevertheless, our simulations pre- 
dict a low but observable level of extended diffuse Ha emis- 
sion in haloes of mass 10 10 /i -1 M©. In fact, extended Ha 
emission forming an Ha-halo, that extends several kilopar- 
secs away from the disc, has been reported for several galax- 
ies. The strength and spatial extent of the emission seem to 
be positive ly correlated with the SFR in the galaxy (e.g., 
iRandll 19961 ) . This highlights the importance of studying CR- 
driven winds including a better description of star formation 
and stellar feedback. 

Some galaxies also exhibit radio haloes with an overall 
correlation between Ha and radio emission. Such a corre- 
lation would naturally arise when CRs are convected into 
the halo, as is the case in our simulations. To discrimi- 
nate a pure advective scenario from our model, the spec- 
tral age of CR ele ctrons could potentially be decisive (see 
iHeesen et al.ll2009l . for such a measurement). Including CR 
streaming should decrease the spectral age since streaming 
electrons escape faster into the halo. However, the adiabatic 
expansion history in both scenarios needs to be identical to 
allow for a fair comparison. Future work that includes MHD 
will address those questions and allow for more definite pre- 
dictions of the expected radio synchrotron emission in this 
scenario. 
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APPENDIX A: IMPLEMENTATION OF CR 
STREAMING 

Here, we discuss our implementation of CR streaming into 
the Lagrangian SPH code Gadget- 2. First, we translate the 
transport equations for CR energy and number density into 
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Lagrangian form and then cast them into the language of 
SPH. Our starting point are the evolution equations for CR 
energy density and number density which can be derive d 
from the CR transport equation (e.g., Skil lindll97lL Il975l ). 



<9£cr 
dt 

dn cr 



(V + Vst) • VPcr - V ' [(V + Vat) (^cr + Per)] , (Al) 
-V- [(V + Vet) Tier], (A2) 



where v is the gas speed, v st is the streaming speed of CRs 
and Per is the CR pressure. Here and in the reminder of this 
derivation, we neglect all CR source terms as well as CR dif- 
fusion. We model CR s ources and sinks explic i tly in our nu- 
merical CR formalism (iPfrommer et al.l 120061 ; lEnfilin et al.l 
20071; ljubelgas et all 120081 ). For the limiting case of strong 
CR-wave scattering (the Bohm limit), the effective diffusion 
velocity is usually much slower than the CR streaming veloc- 
ity, justifying our approximation. As discussed in Section [2 
the streaming speed is assumed to be proportional to the lo- 
cal sound speed, c s , and anti-parallel to a unit vector along 
the CR pressure gradient VP cr in our model, i.e., 



V st 



-Ac s 



VP CI 

IVPci 



(A3) 



where A ^ 1 is a proportionality constant for a plasma where 
the thermal pressure dominates over the magnetic pressure. 
The CR streaming velocity has been taken to be opposite 
to the CR pressure gradient rather than the gradient of the 
CR number density of a CR energy interval (which would 
be the formal criterion for evaluating the direction of the 
CR streaming velocity). For power- law momentum distri- 
butions, ricr oc P cr , and hence the gradients of n cr and P cr 
point into the same direction. If this was not the case, i.e., 
if different energy regimes dominated n cr and P cr , then we 
would be interested in the energy regime that dominates the 
CR pressure (which is responsible for driving the wind) and 
hence in the streaming direction given by the CR pressure 
gradient. 

Since we do not follow MHD in our simulations, we are 
unable to project the CR pressure gradient onto the mag- 
netic field lines to determine the formally correct stream- 
ing direction. However, as we will now argue, our approach 
should correctly capture the main aspects of CR streaming 
and the associated galactic mass loss even without magnetic 
fields. In hydrostatic equilibrium, the gradient of the sum of 
thermal and CR pressure is determined by the gradient of 
the gravitational potential, 

iv(P + P cr ) = -V$. (A4) 

P 

Above the disc, when wave heating induced by streaming 
CRs provides substantial heat input, the CR pressure distri- 
bution determines the total pressure distribution. For CRs 
to escape the galactic disc via streaming requires locally 
open field lines (that are stretched to large heights > 10 kpc 
above the disc). While only 10% of the galactic SN rem- 
nants will create flux tubes that reach those altitudes, SN 
remnant lobes are expected to overlap, suggesting that at 
every location there is suf ficient supply of open field lines 
([Breitschwerdt et al.lll993l ). The field topology at the disc- 
halo interface and beyond will be determined by the buoy- 
ancy of the CR component that drive a larg e-scale u pward 
convection and possibly a galactic mass loss (|Parkerl ll966). 



Hence, the topology of the open field line structure is de- 
termined by the CR pressure stratification which itself is 
shaped by the gravitational potential gradient. This should 
establish a field topology consistent with our assumptions, 
at least when averaged over sufficiently large scales so that 
small-scale fluctuations are averaged out. 

To derive the Lagrangian form of the CR evolution 
equations, we define a Lagrangian time derivative, 

d d „ 

dt = ei+ v ' v > 

and introduce the specific CR energy, s c 
number, n cr , through 






Scrp, 
flcvp. 



(A5) 
and CR particle 

(A6) 

(A7) 



After using the continuity equation for the gas, dp/dt = 
— pV • v, we arrive at the Lagrangian form of the CR evolu- 
tion equations, 

P-^JT = V st • VPcr - PcrV • V - V ■ [v st (pS cr + Per)] ,(A8) 



dncr „ r 

p— 7— = -V • [Vst pn c 



(A9) 



The first term on the right-hand side of equation (|A8[) is the 
wave heating term due to self-excited waves that get (almost 
instantly) damped in the plasma. It results in an energy 
loss of CRs and will be addressed at the end of this Ap- 
pendix. For the moment, we neglect this term in our further 
derivation. The second term on the right-hand side of equa- 
tion (|A8p is the only term involving the plasma velocity v. It 
expresses energy changes due to adiabatic compression (or 
expansion) of CRs when the gas flow converges (diverges). 
Since it is already included in our standard CR implemen- 
tation, we will not discuss it further in the following. The 
remaining term (as well as the only term on the right-hand 
side of equation (|A9|) ) describes CR streaming, i.e., how the 
CR energy and number changes for an individual gas mass 
element due to the CR streaming motion into (or out of) 
that particular mass element. Combining equation (|A3|) for 
the CR streaming velocity with equations (|A8|) and (|A9[) , 
we obtain 

p iif = V -(^ VP ")- ( A1 °) 

p^jL = V-(^VPcr). (All) 

These equations take a form that is similar to a diffusion 
equation with the formal diffusivities 



Ac, 
Ac, 



pScr + Pd 

I VPcr I 
pflcv 

IVPcrl' 



(A12) 
(A13) 



for CR energy density and number density, respectively. 
Therefore, concerning the discretized SPH form of equa- 
tions (IA10[) and (I A 111), the same recipe as f or thermal con- 
ducti on (|Clearvlll999l : ljubelgas et~al1l2004h and CR diffu- 
sion (jJubelgas et al.l 120081 ) can be applieqj. This approach 



6 In Appendix IBll we will show in one-dimensional test cases 
that this indeed leads to reasonable results. 
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rests on an SPH rep resentation of the Laplacian (see, e.g., 
Ijubelgas et al.ll2004r ) that circumvents second derivatives of 
the smoothing kernel, which otherwise would increase nu- 
merical noise. Noise is particularly problematic when an ex- 
plicit time integration scheme is used for the diffusion equa- 
tion, as in our case, because it can cause overshoots that 
may eventually reverse the direction of transport, leading 
to unphysical results unless a quite small timestep is used 
(jjubekas et al.ll2004h . 

Despite using this scheme, we cannot avoid the occur- 
rence of second order kernel derivatives completely, since 
equations (|A12[) and (|A13[) contain the magnitude of the 
CR pressure gradient that also needs to be SPH-smoothed 
and therefore already contains a first-order kernel derivative. 
Leaving this problem aside for the moment, we obtain the 
following evolution equations for the CR quantities of SPH 
particle i, 



OliScr,^ 
~d^ 

dn, 



2£ 



PiPj 



3 ij / J cv,j * cry 



Xij -ViWij, (A14) 



dt 



a PiPj 



where the sum extends over all smoothing neighbours j of 
particle i, and 

,J ,J 

(A16) 



ni + n% 



k% + kI 



(A17) 



represent the harmonic means of the formal diffusiv- 
ities of equations (|A12[) and (|A13[) (each one corre- 
sponding to the positions of i and j), respectively. Fur- 
ther, Wij denotes a symmetrized kernel, i.e., Wij = 
0.5 [W(rij, hi) + W(rij, hj)] and xij — Xi — Xj is the vector 
connecting the positions of both partners. The CR pressure 
gradient at the position of particle i is SPH-smoothed in the 
form 



(•*—! n \ \ ^ / -* cr,i . -* cr,7 

(VP cr ) i = 2^m jPi ^-^- + — Ji 



VWi. 



(A18) 



According to IClearvl ([1999), taking the harmonic mean 
of the formal diffusivities instead of the arithmetic mean 
forces the flux to be continuous in implementations of ther- 
mal conduction, especially in regions where the diffusivity 
changes abruptly. We will see in Appendix [Bl] that this does 
not work perfectly in our case, but irrespectively we apply 
the harmonic mean for the sake of a lower computational 
cost, as will be discussed later on, too. 

If a particle has no CRs at all, its formal diffusivity 
will be zero and therefore also the harmonic means of equa- 
tions (|A12|) and (|A13|) . Thus, there would be no exchange 
between CR pressurized particles and particles without any 
CRs, which is unphysical. In such a case we adopt 

ni = 4, (A19) 

~4 = 4, (A20) 

for the formal diffusivities, i.e., we take the diffusivity of the 
CR pressurized partner. We also employ this choice if the CR 
pressure gradient, VP cr , of one partner is equal to zero, since 
otherwise its formal diffusivity would be infinite. If both 



particles have a CR population and no CR pressure gradient 
in the current time step, then we switch off exchange via 
streaming for this pair. 

Let us now return to the issue of kernel derivatives. 
Indeed, employing equations (|A10[) and (|A11[) for the CR 
streaming implementation will lead to rather noisy results 
because the SPH form of the CR pressure gradient already 
contains one kernel deri vative. To tackl e this problem, we 
follow the suggestions by Ijubelgas et al] (|2004r ) and replace 
the CR pressure of particle j by its SPH-smoothed equiva- 
lent in the evolution equations, i.e., 



E : 



-Lcr.i VVi'i 



(A21) 



where the sum extends over all smoothing neigh bours i of 
partic le j. In agreement with the findings of Ijubelgas et al.l 
(2004), this choice makes our scheme less prone to small- 
scale particle noise, even in the presence of an interpolated 
CR pressure gradient. Thus, our evolution equations for the 
CR quantities now read 

h " - 2 E ^i ( P l~.^ ) *a • V,W«, (A22) 



dt 
drier 



dt 



PiPj 



= ^^'(^fc^) 



Note that we still use the proper CR pressure in the formal 
diffusivity in equation (|A12[h rather than an interpolant. 

Equations (IA23D and (IA22D show that the SPH parti- 
cle i will gain CR energy and CRs from SPH particle j if 
its pressure is less than the (smoothed) pressure of j, con- 
forming with our picture of CRs streaming down their gra- 
dient. Here it is important to note that for particles which 
have no CRs at all, the intrinsic CR pressure, P cr , must not 
be replaced by the smoothed pressure P cr in the transport 
equations above. Otherwise these particles could loose en- 
ergy that they do not have if there are pressurized particles 
in the vicinity that contribute to a non-zero P cr . 

Furthermore, we need to discuss the time step that we 
adopt for our explicit time integration scheme. We decided 
to use a criterion of the form 



1 fm 

At < e— - 

Acs V P 



1/3 



(A24) 



Equation (|A24[) is based on the notion that effects of CR 
streaming should take place on a time-scale equal to the 
time that CRs need to cross the mean inter-particle spacing 
(m/p) ' , travelling at the streaming speed Ac s . Here, e is a 
parameter for which we used e = 0.004 in most of our simu- 
lations. This choice provides us with a good balance between 
accuracy and computational cost. We note that this choice 
would still be fairly expensive for cosmological simulations, 
suggesting that a replacement of the explicit time integra - 
tion by an implicit scheme (e.g., IPetkova fc Springell 12009) 
should be considered in future work. 

In order to ensure conservation of CR energy and par- 
ticle num ber during the excha nge between SPH particles, 
we follow Ijubelgas et al.l (|2008r ) and apply their manifestly 
conservative scheme. Instead of equations (|A23[) and (|A22[) , 
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APPENDIX B: NUMERICAL TESTS 






de c 



dE c 



dt 

drier, 
dt 



dt 
diV cr 



dt 



(A25) 
(A26) 



to calculate the changes of CR energy and particle number 
for an individual SPH particle i due to CR streaming. This 
translates the CR number density and energy density per 
unit mass into proper particle number and energy. Thus, 
the exchange terms between the SPH neighbours i and j 
over which we sum, now read 



dE, 






dt 

dNjj 
dt 



TTliUlj _ij I r C r,j -t cr,i 
PiP 3 S 

TYliTYlj _ij I ±cr,j 
A K,~ 



"ij ' V i W ij , 



(A27) 



Pip j 



Xij-V x Wij. (A28) 



More importantly, in the SPH loop of particle i, we only in- 
ject an amount 0.5 (dEij/dt)AU of energy into its CR pop- 
ulation. In addition, — 0.b(dEij/dt)AU is given to particle 
j at the same time. After j's loop is complete, the effective 
change in the CR energy (and correspondingly in the CR 
particle number) for both particles is then given by 



AEi 



AE, = - 



dE % 



2 [~dt 

fdE 
2 



AU 



dE 



dt J 



dt J 



dE %3 
dt 



AU 



(A29) 
(A30) 



Here we assume for simplicity that j is the only neighbour of 
i and vice versa. In total, during each SPH loop, the change 
in the CR energy and particle number, respectively, is zero 
even if each particle has its own time step. Hence, the total 
change in CR energy and particle number, equations (|A29[) 
and (|A30[) . may be seen as injection with the average time 
step, 0.5 x (Ati + At,), provided that dEij/dt = —dEji/dt. 
So far, we have only addressed the streaming part of 
the CR evolution equations that leads to an exchange of 
CRs between different SPH particles and therefore conserves 
the total energy and particle number in CRs. However, in 
Section [2] we saw that the streaming CRs will excite hydro- 
magnetic waves that get quickly damped in the plasma. This 
gives rise to an additional loss term of CR energy (but not of 
CR particles) for individual SPH particles. We account for 
this wave heating term in each time step At by transferring 
an amount 



r, 



Ac s IVPci 



■At 



(A31) 



of energy (per unit mass) from the CR population of each 
particle i into its thermal reservoir. This will lead to a con- 
tinuous heating of the gas as long the magnitude of the CR 
pressure gradient at the specific location is non-zero, i.e., as 
long as the CRs are streaming through the rest frame of the 
gas. The updated CR spectral cutoff and normalization are 
derived from CR streaming-induced cha nges in Ae and An 
according to the same scheme as in ljubelgas et al.l (|2008r ) . 



Bl Testing the Streaming Implementation 

In this section we report results for our CR streaming imple- 
mentation for a number of one-dimensional test cases. We 
face the difficulty here that in general no analytical solu- 
tions of equations (|A8[) and (|A9[) . respectively, exist. This 
is due to the fact that the advection speed (the streaming 
velocity, equation (|A3|) ) changes sign at extrema in the CR 
pressure distribution. Only in some trivial cases, where no 
extrema are present, it is possible to give an analytical so- 
lution. Therefore, the only way to verify our scheme in a 
more realistic setting is to compare it to results obtained 
with other numerical recipes. There exist several techniques 
to numerically solve advection-type equations with a spa- 
tially varying transport velocity. We decided to use Fromm's 
method which belongs to th e class of f inite v olume methods 
and is described in detail in iLeVequd ([20020 . Also, a recent 
discussion of the properties of the streaming equation as well 
as different ways t o numerically solve it on a one- dimensional 
grid was given bv lSharma et al.l (J2010|). 

As far as the initial conditions for our implementation 
in Gadget- 2 are concerned, we use an ensemble of 10 4 SPH 
particles at rest (i.e. \v\ — 0) in a periodic slab of matter. 
Both, the gravitational and the hydrodynamical accelera- 
tions are switched off to ensure that the changes in the CR 
properties are solely caused by CR streaming. For the same 
reason, we neglect Coulomb, hadronic and wave damping 
losses. The slab has a volume of 10 x 10 x 100 kpc 3 , cor- 
responding to a mean inter-particle spacing of 1 kpc, and 
the particle mass is chosen such that the matter density is 
10 10 M©kpc -3 . Since we use periodic boundary conditions 
for the short dimensions of our simulation box, boundary 
effects are avoided and the setting becomes effectively one- 
dimensional. The temperature throughout the simulation 
box is assumed to be constant and chosen such that the 
sound speed is c s = 20kms _1 . Furthermore, in order to test 
our scheme in the presence of numerical noise in the gas den- 
sity, as encountered in cosmological simulations, the initial 
particle positions resemble an irregular, glass-like distribu- 
tion. The simulations are carried out using 50 smoothing 
neighbours. Regarding the parameters of the CR feedback 
model, we adopt A = 1 and e = 0.004 for the streaming part 
and the CR populations are initialized with a fixed spectral 
slope a = 2.25 and a fixed low-momentum cut-off q = 10.0, 
throughout the simulation volume. Note that the choice for 
e given above leads to a comparatively small time step, 
which is necessary to prevent overshooting in our explicit 
time integration scheme, but results in a high computational 
cost. The expensiveness of numerically solving the stream- 
ing equation was also encountered in the one-dimensional 
tests conducted bv lSharma et al.l (J2010|), even when implicit 
methods were used. The grid for our one-dimensional ref- 
erence simulations adopting Fromm's method is made up 
of 1000 mesh points and are based on a time step of size 
At= 10- 4 /i/^st. 

The streaming equations that we solve read 



de cr 
dt 


1 d 

p dx 


drier 


1 d 


dt 


p dx 



[Vst (p£cr + Per)] , 
(VstpTlcr) , 



(Bl) 

(B2) 
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Figure Bl. Numerical results of the CR streaming implementation in Gadget-2 for the time evolution of an initial Gaussian distribution 
in CR energy (first two rows) and number density (last two rows), respectively. From the top left to the bottom right panel in the first 
two rows, the simulation times increase as t = 0.0, 0.1, 0.2, 0.5, 1.0, 1.5 Gyr, and correspondingly for the last two rows. The black dots 
indicate the SPH particle values, the red line represents the results obtained via Fromm's method, and the dashed line indicates the 
initial distribution. The sound speed is c s = 20kms — 1 throughout the volume, but the propagation speed is larger than that, because of 
the appearance of CR pressure in the formal diffusivities. The SPH particle positions resemble a "glass" -like distribution with a mean 
inter-particle spacing of 1 kpc. For this test, we switched off the gravitational and hydrodynamical accelerations. 



with the magnitude of the streaming speed given by the 
sound speed and the direction determined by the local CR 
pressure gradient, 



^st 



-C s 

+ C S 



for 
for 



dx 
9P„ 



>o, 

<0. 



(B3) 



While the CR description in Gadget- 2 allows for a varying 
effective adiabatic index 7 cr of CRs, we do not account for 
this possibility in our simplified Fromm scheme. Instead, 
we adopt a fixed adiabatic index of 7 cr ~ 1.34, suitable 



for the spectral slope and cut-off given above. Using the 
relation P CY — (7 cr — 1) ps C r, the CR pressure is then simply 
proportional to the energy density and can be replaced by 
the latter in the streaming equation, equation (|B1[) , evolved 
by the Fromm scheme. 

Figure IB1I compares the results obtained with our 
CR streaming implementation and those obtained using 
Fromm's method for the time evolution of the streaming 
equations given above. Here, we simulated the spatial and 
temporal evolution of the CR population, starting from an 
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Figure B2. Resolution study of mass loss history, M(t), in the 10 9 h~ x Mq halo. Left panel: despite varying the number of SPH and 
DM particles used to sample dark matter and baryonic gas by a factor of 20, M(t) is converged. Right panel: we vary the parameter e, 
which regulates the size of the streaming time step of equation (|A24J) , and demonstrate approximate convergence of M(t). 



initial Gaussian distribution in both, CR energy density and 
number density, over a time interval of 1.5 Gyr. 

While the overall agreement of both schemes is good, 
there exist significant differences in the transition region be- 
tween the flat inner part of the CR energy density and the 
outer flanks, with Fromm's method showing a sharper tran- 
sition. This discrepancy is most likely due to the fact that 
SPH is not able to accurately describe sharp discontinu- 
ities but rather tends to smooth them out. Fromm's method 
in contrast should give a numerical solution quite close to 
reality, since the flat inner part of the distribution arises 
owing to the diffusive nature o f the CR streaming e quation 
around extrema, as shown by ISharma et al.l (J2010|) . How- 
ever, in both cases the flanks seem to travel at almost the 
same speed, which indicates a good representation of the 
advective aspect of the CR streaming equations in our im- 
plementation. 

Basically all the aforementioned considerations also 
hold true for the evolution of CR number density, but ad- 
ditional complications arise. Starting at about t — 0.5 Gyr, 
the results become quite noisy. However, this is not problem- 
atic, since the bulk of the SPH particles shows no extreme 
scatter. More importantly, our main requirement is to ob- 
tain a good representation of the dynamically relevant CR 
quantities, i.e., the pressure and energy density. Since this 
is fulfilled, the noise in number density should not signifi- 
cantly influence our results concerning the dynamical impact 
of CRs on the evolution of galactic winds. Note also that the 
propagation of CR number density is slower than that of CR 
energy density. Hence streaming is not able to smooth out 
the former until t = 1.5 Gyr but this is because of the disap- 
pearance of the CR pressure gradient (and correspondingly 
v st = 0) rather than due to inaccuracies in our scheme. 

Apart from an initial Gaussian distribution, we also 
tested a step function with each side of the step set to a 
different sound speed (and thus v s t)- Also in this case we ob- 
served a good overall agreement of CR pressure and energy 
density evolution in both schemes. However, for large steps 
in the sound speed (temperature) the CR energy flux be- 
came discontinuous for several SPH particles located at the 
step. This is not an issue for our isolated haloes but could 



be problematic in simulations of cosmological structure for- 
mation. We have found a slightly modified SPH formulation 
of CR streaming, that behaves better in the aforementioned 
test case. By taking the geometric mean of the formal diffu- 
sivities in the streaming equations, rather than the harmonic 



K % +K? 



\I k % • ni , 



(B4) 



the CR energy flux is well-behaved even at large jumps in 
the temperature field. However, due to the geometric mean 
being bigger than the harmonic one, this modified scheme 
needs a more conservative time step and was therefore not 
used for our isolated halo simulations. In order to reduce 
the computational cost as well as to make the scheme more 
stable for cosmological simulations, we intend to implement 
an implicit time integration scheme in the future. 

B2 Resolution Tests 

We perform a number of resolution tests to study numer- 
ical convergence of the mass loss history and CR Alfven- 
wave heating. In Fig. IB2[ we demonstrate that we obtain 
converged results for our standard choice of the parameter 
e = 0.004, which regulates the size of the streaming time 
step of equation (|A24[) , and the number of SPH and DM 
particles of 1 x 10 5 and 2 x 10 5 , respectively. 
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